##### R CODE TO REPLICATE THE STATISTICAL ANALYSES AND FIGURES OF
##### German Feierherd and Adrian Lucardi, "When the Partisan Becomes Personal: Mayoral Incumbency Effects in Buenos Aires" Forthcoming at the Journal of Elections, Public Opinion and Parties

### The following analysis were run with R version 4.1.2 (2021-11-01) -- "Bird Hippie"
### NOTE: all tables and figures can be searched by name in this file, but they do NOT appear in the order they are included in the text


## emptying the workspace
rm (list=ls ())


## saving the output to a log file
sink (file="Log Analysis When the Partisan Becomes Personal.txt", append=FALSE, type="output", split=FALSE)


## Load/install packages --
if (!require ("pacman")) install.packages ("pacman")
pacman::p_load (
 dplyr, 
 fastDummies,
 ggcorrplot,
 ggplot2,
 ggpmisc,
 ggthemes, ## para theme_map()
 gtools,
 lfe,
 lpdensity,
 readxl,
 rddensity,
 rdrobust,
 schoolmath,
 sf,
 stringr,
 tidyr,
 tidyverse,
 viridis,
 wesanderson,
 xtable )


## function for extracting values of interest from a RD model
extract_rd <- function (x) {
 out <- c (
  as.character (round (x$coef[1], 2))
  , paste0 ("[", round (x$ci[3,1], 2), ":", round (x$ci[3,2], 2), "]")
  , paste0 (x$N_h[1], " | ", x$N_h[2])
  , as.character (round (x$bws[1,1], 1))
  , ""
  )
 return (out) }


## function for extracting values of interest from a RD regression
rd_hte <- function (all, hi, lo, title, party, group1, group2, mean_y_all, mean_y_lo, mean_y_hi){
 
 diff_hilo <- list()
 diff_hilo$est <- hi$coef[1] - lo$coef[1]
 diff_hilo$se <- sqrt((hi$se[3])^2 + (lo$se[3])^2) # 1: conventional; 3: robust 
 diff_hilo$pv <- (1 - pnorm (q=abs(diff_hilo$est / diff_hilo$se)))*2
 depvar <- c (paste0 (title, ",", party), "XXXX", "XXXX", "XXXX")
 model <- c ("ALL", "X=1", "X=0", "DIFF")
 b <- sprintf ("%.2f", c (all$coef[1], hi$coef[1], lo$coef[1], diff_hilo$est[1]))
 ci <- c (paste0 ("[", sprintf ("%.1f", all$ci[2,1]), ":", sprintf ("%.1f", all$ci[2,2]), "]"), 
      paste0 ("[", sprintf ("%.1f", hi$ci[2,1]), ":", sprintf ("%.1f", hi$ci[2,2]), "]"),
      paste0 ("[", sprintf ("%.1f", lo$ci[2,1]), ":", sprintf ("%.1f", lo$ci[2,2]), "]"),
      paste0 ("[", sprintf ("%.1f", diff_hilo$est-qnorm (0.975)*diff_hilo$se), ":",
          sprintf ("%.1f", diff_hilo$est+qnorm (0.975)*diff_hilo$se), "]"))
 pval <- sprintf ("%.2f", c (all$pv[3], hi$pv[3], lo$pv[3], diff_hilo$pv))
 bw <- sprintf ("%.2f", c (all$bws[1,1], hi$bws[1,1], lo$bws[1,1], NA))
 nobs <- c (paste0 (sprintf ("%.0f", all$N_h[1]), "$|$", sprintf ("%.0f", all$N_h[2]))
       , paste0 (sprintf ("%.0f", hi$N_h[1]), "$|$", sprintf ("%.0f", hi$N_h[2]))
       , paste0 (sprintf ("%.0f", lo$N_h[1]), "$|$", sprintf ("%.0f", lo$N_h[2])), NA)
 mean_y <- sprintf ("%.2f", c (mean_y_all, mean_y_hi, mean_y_lo, NA))
 
 tab0 <- rbind (depvar, model, b#, se
         , ci, pval, bw, #ncontrol, ntreat
         nobs, mean_y)
 
 return (tab0)
 
}


## display options
options (digits=4
     , scipen=999 ## disable sci notation
     , show.signif.stars=FALSE
     , max.print=2000
     , tibble.width=Inf
     , tibble.print_max=Inf
     , tibble.print_min=1)


## setting the graphical parameters
(digits <- 2) ## number of decimals to use

col_cutoff <- "black"
col_movie <- "Rushmore1"
wes_palette (col_movie) ## see https://github.com/karthik/wesanderson for a list of palettes
col_bin <- wes_palette (col_movie)[3]
col_line <- wes_palette (col_movie)[5]
alpha_bin <- 1
place_text_x <- -32.5
place_text_y <- 99
size_bin <- 0.85
size_text <- 3.75
theme_set (theme_minimal (base_size=15))

shade_density_ci <- 0.4
shade_density_hist <- 0.25
col_pj <- wes_palette (col_movie)[3]
col_ucr <- wes_palette (col_movie)[4]
col_concurrent <- "gray67"
col_density_left <- wes_palette (col_movie)[3]
col_density_right <- wes_palette (col_movie)[5]
col_density_hist <- wes_palette (col_movie)[4]
alpha_values <- 0.125
map_palette <- "inferno"
jitter_w <- 1/5
jitter_h <- 5
size_tsline <- 1.25


## setting the working director(ies) --> replace this with path to your own working directory
if (grepl ("/Users/gfeierherd", getwd ()) == TRUE){
 home <- "~/Dropbox/Working_Papers/RDD Argentina/code/replication files/" # German's path
 } else {
 home <- "~/Dropbox/Current projects/RDD Argentina/code/replication files/" # Adrian's path
 }
setwd (home)




#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
#### // 1 // DOWNLOADING THE DATA ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#

## if neither the PJ nor the UCR, who wins? 
load ("Elecciones Municipales PBA 1973-2019.RData")

pba_tmp <- pba %>% 
 filter (intendente==1 & midterm == "concurrent" & ucr==0 & pjoficial==0) %>% 
 mutate (
  party = factor (party))

nrow (pba_tmp) ## 154 elections
with (pba_tmp, table (vecinalista))
with (pba_tmp, table (pj))
with (pba_tmp, table (ari))


## main dataset
load ("PBA 1973-2019_full.RData")
pba <- pba %>% mutate (
 small_council = ifelse (concejoSize <= 14, 1, 0)
 
 , winner_t2 = winner_t2*100 
 , winner_t4 = winner_t4*100
 , concejalMaj_t2 = concejalMaj_t2*100
 , concejalMaj_t4 = concejalMaj_t4*100
 
 , winner_lag2 = winner_lag2*100
 , winner_lag4 = winner_lag4*100
 , concejalMaj_lag2 = concejalMaj_lag2*100
 , concejalMaj_lag4 = concejalMaj_lag4*100
 
 , winner_t2_fe = winner_t2_fe*100
 , winner_t4_fe = winner_t4_fe*100
 , concejalMaj_t2_fe = concejalMaj_t2_fe*100
 , concejalMaj_t4_fe = concejalMaj_t4_fe*100
 
 # code to visually separate PJ and UCR in the plot, and to measure vote share in federal elections
 , year_party = ifelse (party=="pj", year - 1/3, year + 1/3)
 , year_presid = ifelse (year %in% c (1983, 1989, seq (1995, 2019, by=4)), 1, 0)
 , fedSh = ifelse (year_presid==1, presiSh, dipuSh) )
summary (pba)


## seccion electoral
seccion <- read_xlsx ("secciones.xlsx") %>% 
 mutate (
  seccion = factor (seccion)
  , conurba = ifelse (seccion %in% c ("1st", "3rd", "8th"), "Conurbano", "Interior")
  , conurba = factor (conurba)
  )
levels (seccion$seccion) ## all's well

# joining
summary (sort (unique (pba$municipio)) %in% sort (unique (seccion$municipio))) ## no discrepancies. Good
pba <- left_join (seccion, pba, by=c ("municipio" = "municipio")) %>% 
 mutate (
  seccion_party = ifelse (party=="pj", as.numeric (seccion) - 1/8, as.numeric (seccion) + 1/8)
  , year_in2 = year + 2
  , year_in4 = year + 4 ## we'll use these to measure presidential copartisanship + approval in t+2 and t+4
  )


## copartisanship
copartisan <- read_xlsx ("copartisans.xlsx", sheet="Sheet1")

# joining
pba <- left_join (pba, copartisan, by="year") %>% 
 mutate (
  copartisan_presi = ifelse (party == party_presi, 1, 0)
  , copartisan_presi_t2 = ifelse (party == party_presi_t2, 1, 0)
  , copartisan_presi_t4 = ifelse (party == party_presi_t4, 1, 0)
  , copartisan_gover = ifelse (party == party_gover, 1, 0)
  , copartisan_gover_t2 = ifelse (party == party_gover_t2, 1, 0)
  , copartisan_gover_t4 = ifelse (party == party_gover_t4, 1, 0) ) %>% 
  group_by (municipio, party) %>% 
  mutate (
    copartisan_presi_lag2 = lag (copartisan_presi, n=1, order_by=year)
    , copartisan_presi_lag4 = lag (copartisan_presi, n=2, order_by=year)
    )


## executive approval, by quarter
ead <- read_csv ("EAD+2.0+quarter+101019.csv") %>% 
 filter (
  Country == "Argentina") %>% 
 mutate (
  approval = NET_Not_Smoothed
  , election = ifelse (
   (year %in% c (1983, 1985, 1993, 1997, 1999, 2001, 2005, 2007, 2011, 2013, 2015, 2017, 2019) & quarter == 3)
   | (year %in% c (1987, 1991, 2003) & quarter == 2)
   | (year %in% c (1989, 1995, 2009) & quarter == 1)
   , 1, 0) )
median (ead$approval, na.rm=TRUE) ## 4.48
median (filter (ead, election==1)$approval, na.rm=TRUE) ## 3.93. We'll use the latter, though it makes little difference in practice
ead <- ead %>% mutate (
 approval_high = ifelse (approval > median (ead$approval, na.rm=TRUE), 1, 0))
ead %>% filter (election==1) %>% select (year, quarter, qtr, approval_high, approval, Approval_Not_Smoothed)

# there is a large difference, not only in means (-5.3 vs 24.0), but also between the highest value of one case and the lowest of the other (3.9 vs 11.0):
ead %>% filter (election==1) %>% group_by (approval_high) %>% summarise (
 approval_mean = mean (approval, na.rm=TRUE)
 , approval_min = min (approval, na.rm=TRUE)
 , approval_max = max (approval, na.rm=TRUE))

# joining
pba <- left_join (pba
          , select (filter (ead, election==1), year, approval, approval_high)
          , by=c ("year" = "year"))
pba <- left_join (pba
          , select (filter (ead, election==1), year, approval, approval_high)
          , by=c ("year_in2" = "year"), suffix=c ("", "_t2"))
pba <- left_join (pba
                  , select (filter (ead, election==1), year, approval, approval_high)
                  , by=c ("year_in4" = "year"), suffix=c ("", "_t4"))
pba <- pba %>% 
  group_by (municipio, party) %>% 
  mutate (
    approval_high_lag2 = lag (approval_high, n=1, order_by=year)
    , approval_high_lag4 = lag (approval_high, n=2, order_by=year)
  ) %>% 
 mutate (
  favored_presi = ifelse (is.na (approval_high), NA, ifelse (copartisan_presi==1 & approval_high==1, 1, 0))
  , favored_presi_t2 = ifelse (is.na (approval_high_t2), NA, ifelse (copartisan_presi_t2==1 & approval_high_t2==1, 1, 0))
  , favored_presi_t4 = ifelse (is.na (approval_high_t4), NA, ifelse (copartisan_presi_t4==1 & approval_high_t4==1, 1, 0))
  , disfavored_presi = ifelse (is.na (approval_high), NA, ifelse (copartisan_presi==1 & approval_high==0, 1, 0))
  , disfavored_presi_t2 = ifelse (is.na (approval_high_t2), NA, ifelse (copartisan_presi_t2==1 & approval_high_t2==0, 1, 0))
  , disfavored_presi_t4 = ifelse (is.na (approval_high_t4), NA, ifelse (copartisan_presi_t4==1 & approval_high_t4==0, 1, 0))
  , favored_opp = ifelse (is.na (approval_high), NA, ifelse (copartisan_presi==0 & approval_high==0, 1, 0))
  , favored_opp_t2 = ifelse (is.na (approval_high_t2), NA, ifelse (copartisan_presi_t2==0 & approval_high_t2==0, 1, 0))
  , favored_opp_t4 = ifelse (is.na (approval_high_t4), NA, ifelse (copartisan_presi_t4==0 & approval_high_t4==0, 1, 0))
  , favored_presi_lag2 = ifelse (is.na (approval_high_lag2), NA, ifelse (copartisan_presi_lag2==1 & approval_high_lag2==1, 1, 0))
  , favored_presi_lag4 = ifelse (is.na (approval_high_lag4), NA, ifelse (copartisan_presi_lag4==1 & approval_high_lag4==1, 1, 0))
  , disfavored_presi_lag2 = ifelse (is.na (approval_high_lag2), NA, ifelse (copartisan_presi_lag2==1 & approval_high_lag2==0, 1, 0))
  , disfavored_presi_lag4 = ifelse (is.na (approval_high_lag4), NA, ifelse (copartisan_presi_lag4==1 & approval_high_lag4==0, 1, 0))
  )
summary (pba)
with (pba, table (year, approval_high_t2, useNA="always")) ## perfect


## splitting between PJ and UCR
dpj <- filter (pba, year >= 1983 & party == "pj") ## the "toptwo" variable doesn't matter b/c running is defined as NA if toptwo==0
ducr <- filter (pba, year >= 1983 & party == "ucr")

# in both cases, (dis)favorability status is highly correlated with (a) decade; and (b) whether an election is a midterm one:
dpj %>% filter (midterm=="concurrent" & year <= 2015) %>% select (year, copartisan_presi:copartisan_presi_t4, favored_presi:favored_opp_t4) %>% unique ()
ducr %>% filter (midterm=="concurrent" & year <= 2015) %>% select (year, copartisan_presi:copartisan_presi_t4, favored_presi:favored_opp_t4) %>% unique ()



### some descriptive stats ####

# number of mayoral elections
nrow (select (filter (pba, intendente==1), municipio, year)) ## 1175

# number of registered voters in median municipality (2019)
quantile (filter (pba, year==2019 & party=='pj')$padron, probs=c (0.05, 0.25, 0.5, 2/3, 0.75, 0.95))
quantile (filter (pba, year==2019 & party=='ucr')$padron, probs=c (0.05, 0.25, 0.5, 2/3, 0.75, 0.95))
filter (pba, year==2019 & party=='pj' & municipio=="LA MATANZA")$padron ## 1.117M
filter (pba, year==2019 & party=='pj' & municipio=="LOMAS DE ZAMORA")$padron ## 543K

# mayoral re-runs and re-elections (2011 & 2015 cohorts): Table A2
tab1 <- read_excel ("PBA_2011_2019.xlsx")
summary (tab1)

d1 <- tab1 %>% 
  group_by (year) %>% 
  select (year, reruns, reelects) %>% 
  summarize(
    n_munis = n ()
    , count_rerun = sum (reruns)
    , share_rerun = 100*mean (reruns)
    , count_reelects = sum (reelects)
    , share_reelects = 100*count_reelects/count_rerun
  )

d2 <- tab1 %>% 
  select(reruns,reelects) %>% 
  summarize(
    n_munis = n ()
    , count_rerun = sum(reruns)
    , share_rerun = 100*mean(reruns)
    , count_reelects = sum(reelects)
    , share_reelects = 100*count_reelects/count_rerun
  ) %>% 
  mutate(year="Total")

(d12 = rbind(d1, d2) )
names(d12) <- c("", "sample size (N)", "Re-run (N)", "Re-run (\\%)", "Re-elects (N)", "Re-elects (\\%)") 

print (xtable (d12, type = "latex", align=c("l","c","c","c","c","c","c"), 
               caption="Re-running and re-election rates among PBA mayors (2011, 2015)", 
               label="T:re"), include.rownames=FALSE)


## average turnout (1993-2019 only)
summary (filter (pba, party=="pj" & !is.na (padron))$year) ## 1993 to 2019
with (filter (pba, party=="pj" & !is.na (padron)), mean (emitidos/padron)*100) ## 80.82%
with (filter (pba, party=="pj" & !is.na (padron)), by (emitidos/padron*100, midterm, mean)) ## 2pp higher in concurrent elections (81.88% vs. 79.5%)


## % of elections in which the PJ and UCR participated (and won)

# Between 1983 and 2019, there was a total of 1319 elections in concurrent years:
(tot_ele <- nrow (unique (select (filter (pba, year %in% seq (1983, 2019, by=4)), municipio, year))))
(tot_pj <- nrow (unique (select (filter (pba, party=='pj' & year %in% seq (1983, 2019, by=4)), municipio, year)))); round (tot_pj/tot_ele*100, 2) ## PJ participated in 1314
(tot_ucr <- nrow (unique (select (filter (pba, party=='ucr' & year %in% seq (1983, 2019, by=4)), municipio, year)))); round (tot_ucr/tot_ele*100, 2) ## UCR participated in 1274

(win_both <- sum (filter (pba, year %in% seq (1983, 2019, by=4))$winner)); round (win_both/tot_ele*100, 2) ## won 88.25 between the two





#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
#### // 2 // RD PLOTS, MAIN ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#

## updating WD to export the figures there
setwd (str_c (home, "figures/"))


## list of outcome variables + objects to store the results
out_main <- c ("winner_t2", "winner_t4", "voteSh_t2", "voteSh_t4", "concejalSh_t2", "concejalSh_t4", "dipuSh_t2", "presiSh_t4")
out_lag <- c ("winner_lag2", "winner_lag4", "voteSh_lag2", "voteSh_lag4", "concejalSh_lag2", "concejalSh_lag4", "dipuSh_lag2", "presiSh_lag4")
out_fe <- c ("winner_t2_fe", "winner_t4_fe", "voteSh_t2_fe", "voteSh_t4_fe", "concejalSh_t2_fe", "concejalSh_t4_fe", "dipuSh_t2_fe", "presiSh_t4_fe")
out_midt <- c ("winner_t2", "winner_t4", "voteSh_t2", "voteSh_t4", "concejalSh_t2", "concejalSh_t4", "presiSh_t2", "dipuSh_t4")
outcomes <- c (out_main, out_lag, out_fe)
plot_results <- NULL ## we need this to "store" the results


### (2.1) loop to get the underlying data ####
for (o in 1:length (outcomes)){
 
 ### selecting the outcomes of interest
 dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == outcomes[o])])
 ducr$outcome_tmp <- unlist (ducr[,which (colnames (ducr) == outcomes[o])])
 
 
 ### drawing the plots
 
 ## (2.1.1) running measured in concurrent election (i.e., Treatment is mayoral incumbency)
 
 # QSMV: quantile-spaced, mimicking variability -> better reflects the actual distribution of the data
 qsmv_u3_pj11 <- with (filter (dpj, midterm=="concurrent"), rdplot (
  y=outcome_tmp, x=running, p=3, binselect="qsmv", scale=1, kernel="uniform", x.lim=c (-50, 50)))
 qsmv_u3_ucr11 <- with (filter (ducr, midterm=="concurrent"), rdplot (
  y=outcome_tmp, x=running, p=3, binselect="qsmv", scale=1, kernel="uniform", x.lim=c (-50, 50)))
 
 # exporting the data we'll use to build the plots manually later
 qsmv_u3_pj11p <- bind_rows (
  ggplot_build (qsmv_u3_pj11$rdplot)$data[[1]] %>% mutate (type="bin")
  , ggplot_build (qsmv_u3_pj11$rdplot)$data[[2]] %>% mutate (type="ll")
  , ggplot_build (qsmv_u3_pj11$rdplot)$data[[3]] %>% mutate (type="rl") ) %>% mutate (
   bin_select="qsmv", poly=3, party="pj", midterm="concurrent", sample="full", outcome=outcomes[o] )
 
  qsmv_u3_ucr11p <- bind_rows (
  ggplot_build (qsmv_u3_ucr11$rdplot)$data[[1]] %>% mutate (type="bin")
  , ggplot_build (qsmv_u3_ucr11$rdplot)$data[[2]] %>% mutate (type="ll")
  , ggplot_build (qsmv_u3_ucr11$rdplot)$data[[3]] %>% mutate (type="rl") ) %>% mutate (
   bin_select="qsmv", poly=3, party="ucr", midterm="concurrent", sample="full", outcome=outcomes[o] )
 
 
 ## (2.1.2) running measured in MIDTERM election (i.e., Treatment is being the most voted party in a midterm, which grants no special status)
 
 # QSMV: quantile-spaced, mimicking variability
 qsmv_u3_pj21 <- with (filter (dpj, midterm=="midterm"), rdplot (
  y=outcome_tmp, x=running, p=3, binselect="qsmv", scale=1, kernel="uniform", x.lim=c (-50, 50)))
 qsmv_u3_ucr21 <- with (filter (ducr, midterm=="midterm"), rdplot (
  y=outcome_tmp, x=running, p=3, binselect="qsmv", scale=1, kernel="uniform", x.lim=c (-50, 50)))
 
 # exporting the data we'll use to build the plots manually later
 qsmv_u3_pj21p <- bind_rows (
  ggplot_build (qsmv_u3_pj21$rdplot)$data[[1]] %>% mutate (type="bin")
  , ggplot_build (qsmv_u3_pj21$rdplot)$data[[2]] %>% mutate (type="ll")
  , ggplot_build (qsmv_u3_pj21$rdplot)$data[[3]] %>% mutate (type="rl") ) %>% mutate (
   bin_select="qsmv", poly=3, party="pj", midterm="midterm", sample="full", outcome=outcomes[o] )
 qsmv_u3_ucr21p <- bind_rows (
  ggplot_build (qsmv_u3_ucr21$rdplot)$data[[1]] %>% mutate (type="bin")
  , ggplot_build (qsmv_u3_ucr21$rdplot)$data[[2]] %>% mutate (type="ll")
  , ggplot_build (qsmv_u3_ucr21$rdplot)$data[[3]] %>% mutate (type="rl") ) %>% mutate (
   bin_select="qsmv", poly=3, party="ucr", midterm="midterm", sample="full", outcome=outcomes[o] )
 
 
 ## exporting the results: no controls
 plot_results <- bind_rows (plot_results,
             qsmv_u3_pj11p, qsmv_u3_ucr11p
             , qsmv_u3_pj21p, qsmv_u3_ucr21p) }


## updating some variables
plot_results2 <- plot_results %>%
 separate (outcome, c ("out_var", "out_time", "out_fe"), remove=F, sep="_") %>% 
  mutate (
 type = factor (type)
 , bin_select = factor (bin_select)
 , party = factor (party)
 , midterm = factor (midterm)
 , sample = factor (sample)
 , outcome = factor (outcome)
 , out_type = ifelse (outcome %in% out_main, "main", ifelse (
  outcome %in% out_lag, "placebo - lagged", "main - FE"))
 , out_type = factor (out_type, levels = c ("main", "main - FE", "placebo - lagged"))
 
 # out_var: type of outcome variable (victory, vote share, etc)
 , out_var = recode (out_var
           , winner = "victory (0/100)"
           , voteSh = "vote share\n(municipal) (0:100)"
           , concejalSh = "seat share\n(municipal) (0:100)"
           , concejalMaj = "majority in\ncouncil (0/100)"
           , presiSh = "vote share\n(federal) (0:100)"
           , dipuSh = "vote share\n(federal) (0:100)" )
 , out_var = factor (out_var, levels = c (
  "victory (0/100)", "vote share\n(municipal) (0:100)", "seat share\n(municipal) (0:100)", "majority in\ncouncil (0/100)", "vote share\n(federal) (0:100)"))
 
 # out_var: time at which outcome variable is measured (t+2, t+4, t-2, etc)
 , out_time = recode (out_time, lag4 = "t - 4", lag2 = "t - 2"
            , t2 = "t + 2", t4 = "t + 4")
 , out_time = factor (out_time)
 , out_fe = factor (out_fe) ) %>% group_by (
  party, midterm, sample, outcome, bin_select, poly) %>% mutate ( ## to identify each plot uniquely
   y_below = ifelse (type=="ll" & x==0, y, NA)
   , y_below = max (y_below, na.rm=T)
   , y_above = ifelse (type=="rl" & x==0, y, NA)
   , y_above = max (y_above, na.rm=T)
   , rd_estim = y_above - y_below
   , rd_estim_text = paste0 ("\u03b2 = ", sprintf ("%.1f", round (rd_estim, 1)), "pp.")
   
   , PANEL = NULL
   , group = NULL
   , shape = NULL
   , colour = NULL
   , size = NULL
   , fill = NULL
   , alpha = NULL
   , stroke = NULL
   , flipped_aes = NULL
   , linetype = NULL ) ## warning reported. Don't worry
summary (plot_results2)



### (2.2) building the plots proper ####

## (2.2.1) Main results: effect of incumbency on 8 outcome variables ####

# (2.2.1.1) PJ (Figure 1a)
(rdp11_mainpj_qsmv3 <- ggplot (filter (plot_results2, out_type=='main', sample=="full", midterm=="concurrent", party=="pj", bin_select=="qsmv", poly==3, type=="bin"), aes (x=x, y=y))
 + geom_vline (xintercept=0, col=col_cutoff)
 + geom_point (size=size_bin, alpha=alpha_bin, col=col_bin)
 + geom_line (data=filter (plot_results2, out_type=='main', sample=="full", midterm=="concurrent", party=="pj", type=="ll"), aes (x=x, y=y), col=col_line)
 + geom_line (data=filter (plot_results2, out_type=='main', sample=="full", midterm=="concurrent", party=="pj", type=="rl"), aes (x=x, y=y), col=col_line)
 + facet_grid (out_time ~ out_var)
 + geom_text (data=filter (plot_results2, out_type=='main', sample=="full", midterm=="concurrent", party=="pj", bin_select=="qsmv", poly==3) %>% group_by (out_var, out_time) %>% summarise (rd_estim_text=unique (rd_estim_text))
        , mapping=aes (x=place_text_x, y=place_text_y, label=rd_estim_text), size=size_text )
 + xlim (-50, 50) + ylim (0, 100)
 + xlab (expression (margin~of~victory[t]~(concurrent~election))) + ylab ("")
 + theme (axis.title.x = element_text (margin=margin (t=10, r=0, b=0, l=0))) )

# (2.2.1.2) UCR (Figure 1b)
(rdp21_mainucr_qsmv3 <- ggplot (filter (plot_results2, out_type=='main', sample=="full", midterm=="concurrent", party=="ucr", bin_select=="qsmv", poly==3, type=="bin"), aes (x=x, y=y))
 + geom_vline (xintercept=0, col=col_cutoff)
 + geom_point (size=size_bin, alpha=alpha_bin, col=col_bin)
 + geom_line (data=filter (plot_results2, out_type=='main', sample=="full", midterm=="concurrent", party=="ucr", type=="ll"), aes (x=x, y=y), col=col_line)
 + geom_line (data=filter (plot_results2, out_type=='main', sample=="full", midterm=="concurrent", party=="ucr", type=="rl"), aes (x=x, y=y), col=col_line)
 + facet_grid (out_time ~ out_var)
 + geom_text (data=filter (plot_results2, out_type=='main', sample=="full", midterm=="concurrent", party=="ucr", bin_select=="qsmv", poly==3) %>% group_by (out_var, out_time) %>% summarise (rd_estim_text=unique (rd_estim_text))
        , mapping=aes (x=place_text_x, y=place_text_y, label=rd_estim_text), size=size_text )
 + xlim (-50, 50) + ylim (0, 100)
 + xlab (expression (margin~of~victory[t]~(concurrent~election))) + ylab ("")
 + theme (axis.title.x = element_text (margin=margin (t=10, r=0, b=0, l=0))) )


## (2.2.2) Main results + FE: outcome demeaned of municipality- and year- effects ####

# range of the y axis
summary (filter (plot_results2, out_type=="main - FE")$y) ## -37 to 55; we make it -40 to 60. We also change the location of the text

# (2.2.2.1) PJ (Figure A6a)
(rdp12_fepj_qsmv3 <- ggplot (filter (plot_results2, out_type=='main - FE', sample=="full", midterm=="concurrent", party=="pj", bin_select=="qsmv", poly==3, type=="bin"), aes (x=x, y=y))
 + geom_vline (xintercept=0, col=col_cutoff)
 + geom_point (size=size_bin, alpha=alpha_bin, col=col_bin)
 + geom_line (data=filter (plot_results2, out_type=='main - FE', sample=="full", midterm=="concurrent", party=="pj", type=="ll"), aes (x=x, y=y), col=col_line)
 + geom_line (data=filter (plot_results2, out_type=='main - FE', sample=="full", midterm=="concurrent", party=="pj", type=="rl"), aes (x=x, y=y), col=col_line)
 + facet_grid (out_time ~ out_var)
 + geom_text (data=filter (plot_results2, out_type=='main - FE', sample=="full", midterm=="concurrent", party=="pj", bin_select=="qsmv", poly==3) %>% group_by (out_var, out_time) %>% summarise (rd_estim_text=unique (rd_estim_text))
        , mapping=aes (x=place_text_x, y=place_text_y-60, label=rd_estim_text), size=size_text )
 + xlim (-50, 50) + ylim (-40, 40)
 + xlab (expression (margin~of~victory[t]~(concurrent~election))) + ylab ("")
 + theme (axis.title.x = element_text (margin=margin (t=10, r=0, b=0, l=0))) )

# (2.2.2.2) UCR (Figure A6b)
(rdp22_feucr_qsmv3 <- ggplot (filter (plot_results2, out_type=='main - FE', sample=="full", midterm=="concurrent", party=="ucr", bin_select=="qsmv", poly==3, type=="bin"), aes (x=x, y=y))
 + geom_vline (xintercept=0, col=col_cutoff)
 + geom_point (size=size_bin, alpha=alpha_bin, col=col_bin)
 + geom_line (data=filter (plot_results2, out_type=='main - FE', sample=="full", midterm=="concurrent", party=="ucr", type=="ll"), aes (x=x, y=y), col=col_line)
 + geom_line (data=filter (plot_results2, out_type=='main - FE', sample=="full", midterm=="concurrent", party=="ucr", type=="rl"), aes (x=x, y=y), col=col_line)
 + facet_grid (out_time ~ out_var)
 + geom_text (data=filter (plot_results2, out_type=='main - FE', sample=="full", midterm=="concurrent", party=="ucr", bin_select=="qsmv", poly==3) %>% group_by (out_var, out_time) %>% summarise (rd_estim_text=unique (rd_estim_text))
        , mapping=aes (x=place_text_x, y=place_text_y-60, label=rd_estim_text), size=size_text )
 + xlim (-50, 50) + ylim (-40, 40)
 + xlab (expression (margin~of~victory[t]~(concurrent~election))) + ylab ("")
 + theme (axis.title.x = element_text (margin=margin (t=10, r=0, b=0, l=0))) )


## (2.2.3) Placebo (I): there should be no effect on lagged DVs ####

# (2.2.3.1) PJ (Figure A9a)
(rdp13_lagpj_qsmv3 <- ggplot (filter (plot_results2, out_type=='placebo - lagged', sample=="full", midterm=="concurrent", party=="pj", bin_select=="qsmv", poly==3, type=="bin"), aes (x=x, y=y))
 + geom_vline (xintercept=0, col=col_cutoff)
 + geom_point (size=size_bin, alpha=alpha_bin, col=col_bin)
 + geom_line (data=filter (plot_results2, out_type=='placebo - lagged', sample=="full", midterm=="concurrent", party=="pj", type=="ll"), aes (x=x, y=y), col=col_line)
 + geom_line (data=filter (plot_results2, out_type=='placebo - lagged', sample=="full", midterm=="concurrent", party=="pj", type=="rl"), aes (x=x, y=y), col=col_line)
 + facet_grid (out_time ~ out_var)
 + geom_text (data=filter (plot_results2, out_type=='placebo - lagged', sample=="full", midterm=="concurrent", party=="pj", bin_select=="qsmv", poly==3) %>% group_by (out_var, out_time) %>% summarise (rd_estim_text=unique (rd_estim_text))
        , mapping=aes (x=place_text_x, y=place_text_y, label=rd_estim_text), size=size_text )
 + xlim (-50, 50) + ylim (0, 100)
 + xlab (expression (margin~of~victory[t]~(concurrent~election))) + ylab ("")
 + theme (axis.title.x = element_text (margin=margin (t=10, r=0, b=0, l=0))) )

# (2.2.3.2) UCR (Figure A9b)
(rdp23_lagucr_qsmv3 <- ggplot (filter (plot_results2, out_type=='placebo - lagged', sample=="full", midterm=="concurrent", party=="ucr", bin_select=="qsmv", poly==3, type=="bin"), aes (x=x, y=y))
 + geom_vline (xintercept=0, col=col_cutoff)
 + geom_point (size=size_bin, alpha=alpha_bin, col=col_bin)
 + geom_line (data=filter (plot_results2, out_type=='placebo - lagged', sample=="full", midterm=="concurrent", party=="ucr", type=="ll"), aes (x=x, y=y), col=col_line)
 + geom_line (data=filter (plot_results2, out_type=='placebo - lagged', sample=="full", midterm=="concurrent", party=="ucr", type=="rl"), aes (x=x, y=y), col=col_line)
 + facet_grid (out_time ~ out_var)
 + geom_text (data=filter (plot_results2, out_type=='placebo - lagged', sample=="full", midterm=="concurrent", party=="ucr", bin_select=="qsmv", poly==3) %>% group_by (out_var, out_time) %>% summarise (rd_estim_text=unique (rd_estim_text))
        , mapping=aes (x=place_text_x, y=place_text_y, label=rd_estim_text), size=size_text )
 + xlim (-50, 50) + ylim (0, 100)
 + xlab (expression (margin~of~victory[t]~(concurrent~election))) + ylab ("")
 + theme (axis.title.x = element_text (margin=margin (t=10, r=0, b=0, l=0))) )


## (2.2.4) Placebo (II): being the most voted party at midterm, which confers no special status, should not matter much ####

# (2.2.4.1) PJ (Figure A7a)
(rdp14_midtpj_qsmv3 <- ggplot (filter (plot_results2, out_type=='main', sample=="full", midterm=="midterm", party=="pj", bin_select=="qsmv", poly==3, type=="bin"), aes (x=x, y=y))
 + geom_vline (xintercept=0, col=col_cutoff)
 + geom_point (size=size_bin, alpha=alpha_bin, col=col_bin)
 + geom_line (data=filter (plot_results2, out_type=='main', sample=="full", midterm=="midterm", party=="pj", type=="ll"), aes (x=x, y=y), col=col_line)
 + geom_line (data=filter (plot_results2, out_type=='main', sample=="full", midterm=="midterm", party=="pj", type=="rl"), aes (x=x, y=y), col=col_line)
 + facet_grid (out_time ~ out_var)
 + geom_text (data=filter (plot_results2, out_type=='main', sample=="full", midterm=="midterm", party=="pj", bin_select=="qsmv", poly==3) %>% group_by (out_var, out_time) %>% summarise (rd_estim_text=unique (rd_estim_text))
        , mapping=aes (x=place_text_x, y=place_text_y, label=rd_estim_text), size=size_text )
 + xlim (-50, 50) + ylim (0, 100)
 + xlab (expression (margin~of~victory[t]~(midterm~election))) + ylab ("")
 + theme (axis.title.x = element_text (margin=margin (t=10, r=0, b=0, l=0))) )

# (2.2.4.2) UCR (Figure A7b)
(rdp24_midtucr_qsmv3 <- ggplot (filter (plot_results2, out_type=='main', sample=="full", midterm=="midterm", party=="ucr", bin_select=="qsmv", poly==3, type=="bin"), aes (x=x, y=y))
 + geom_vline (xintercept=0, col=col_cutoff)
 + geom_point (size=size_bin, alpha=alpha_bin, col=col_bin)
 + geom_line (data=filter (plot_results2, out_type=='main', sample=="full", midterm=="midterm", party=="ucr", type=="ll"), aes (x=x, y=y), col=col_line)
 + geom_line (data=filter (plot_results2, out_type=='main', sample=="full", midterm=="midterm", party=="ucr", type=="rl"), aes (x=x, y=y), col=col_line)
 + facet_grid (out_time ~ out_var)
 + geom_text (data=filter (plot_results2, out_type=='main', sample=="full", midterm=="midterm", party=="ucr", bin_select=="qsmv", poly==3) %>% group_by (out_var, out_time) %>% summarise (rd_estim_text=unique (rd_estim_text))
        , mapping=aes (x=place_text_x, y=place_text_y, label=rd_estim_text), size=size_text )
 + xlim (-50, 50) + ylim (0, 100)
 + xlab (expression (margin~of~victory[t]~(midterm~election))) + ylab ("")
 + theme (axis.title.x = element_text (margin=margin (t=10, r=0, b=0, l=0))) )



### (2.3) exporting the plots ####

## (2.3.1) larger plots (2 per page)

pwid <- 17*1.5
phei <- 17*0.825

# full sample
ggsave ("figRDPlotFullPJ.png"
    , rdp11_mainpj_qsmv3, width=pwid, height=phei, units="cm")
ggsave ("figRDPlotFullUCR.png"
    , rdp21_mainucr_qsmv3, width=pwid, height=phei, units="cm")

# demeaned outcomes ('FEs')
ggsave ("figRDPlotFEPJ.png"
    , rdp12_fepj_qsmv3, width=pwid, height=phei, units="cm")
ggsave ("figRDPlotFEUCR.png"
    , rdp22_feucr_qsmv3, width=pwid, height=phei, units="cm")

# placebo: most voted at midterm
ggsave ("figRDPlotMidtPJ.png"
    , rdp14_midtpj_qsmv3, width=pwid, height=phei, units="cm")
ggsave ("figRDPlotMidtUCR.png"
    , rdp24_midtucr_qsmv3, width=pwid, height=phei, units="cm")

# lagged outcomes
ggsave ("figRDPlotFullLagPJ.png"
    , rdp13_lagpj_qsmv3, width=pwid, height=phei, units="cm")
ggsave ("figRDPlotFullLagUCR.png"
    , rdp23_lagucr_qsmv3, width=pwid, height=phei, units="cm")





#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
#### // 3 // MAIN ANALYSIS ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#

### (3.1) Main results: full sample; and disaggregating by (i) copartisan president; and (ii) whether the copartisan president is favored in the polls (PJ only)

## doing the analysis
rd_results_copartisan <- NULL; rd_results_favored <- NULL## we need this to "store" the results
for (i in 1:length (out_main)){
  
  ## PJ
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_main[i])])
  
  all <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_t2", out_main)) { ## outcome measured at t+2
    hi <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  }
  
  pj_copartisan <- rd_hte (all, hi, lo, out_main[i], "PJ", "Copartisan", "Opposition"
                           , mean_y_all = mean (filter (dpj, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                           , mean_y_hi = ifelse (
                             i %in% grep ("_t2", out_main)
                             , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                             , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                           , mean_y_lo = ifelse (
                             i %in% grep ("_t2", out_main)
                             , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                             , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## UCR
  ducr$outcome_tmp <- unlist (ducr[,which (colnames (ducr) == out_main[i])])
  
  all <- with (filter (ducr, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_t2", out_main)) { ## outcome measured at t+2
    hi <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  }
  
  ucr_copartisan <- rd_hte (all, hi, lo, out_main[i], "UCR", "Copartisan", "Opposition"
                            , mean_y_all = mean (filter (ducr, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                            , mean_y_hi = ifelse (
                              i %in% grep ("_t2", out_main)
                              , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                              , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                            , mean_y_lo = ifelse (
                              i %in% grep ("_t2", out_main)
                              , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                              , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_copartisan <- rbind (rd_results_copartisan, cbind (matrix (unlist (
    list (pj_copartisan, ucr_copartisan)), nrow=8)))
}
rd_results_copartisan


## favored copartisan president: PJ only
for (i in 1:length (out_main)){
  
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_main[i])])
  
  all <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_t2", out_main)) { ## outcome measured at t+2
    hi <- with (filter (dpj, midterm=="concurrent" & favored_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & disfavored_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (dpj, midterm=="concurrent" & favored_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & disfavored_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  }
  
  pj_favored <- rd_hte (all, hi, lo, out_main[i], "PJ", "Favored", "Disfavored"
                        , mean_y_all = mean (filter (dpj, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                        , mean_y_hi = ifelse (
                          i %in% grep ("_t2", out_main)
                          , mean (filter (dpj, midterm=="concurrent" & favored_presi_t2==1 & !is.na (favored_presi_t2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                          , mean (filter (dpj, midterm=="concurrent" & favored_presi_t4==1 & !is.na (favored_presi_t4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                        , mean_y_lo = ifelse (
                          i %in% grep ("_t2", out_main)
                          , mean (filter (dpj, midterm=="concurrent" & disfavored_presi_t2==0 & !is.na (disfavored_presi_t2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                          , mean (filter (dpj, midterm=="concurrent" & disfavored_presi_t4==0 & !is.na (disfavored_presi_t4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_favored <- rbind (rd_results_favored, cbind (matrix (unlist (
    list (pj_favored)), nrow=8)))
}
rd_results_favored


## (3.1.1) full sample: constructing and exporting Table 1 ####
tabMain <- rbind (
 matrix (rd_results_copartisan[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),1]
         , byrow=F, nrow=6)
 , matrix (rd_results_copartisan[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),5]
           , byrow=F, nrow=6))
(tabMain <- cbind (tabMain[,1:2], rep ("", nrow (tabMain)), tabMain[,3:4], rep ("", nrow (tabMain)), tabMain[,5:6], rep ("", nrow (tabMain)), tabMain[,7:8]))

rows <- rep ( c ("estimate ($\\hat{\\tau}_{\\textsc{rd}}$)"
         , "95\\% \\textsc{ci}"
         , "$p$-value"
         , "bwd."
         , "$N$"
         , "control mean"), 2)

Header1 <- paste ("\\toprule & \\multicolumn{2}{c}{} & & \\multicolumn{2}{c}{\\emph{vote share}} & & \\multicolumn{2}{c}{\\emph{seat share}} & & \\multicolumn{2}{c}{\\emph{vote share}} \\\\ \n")
Header2 <- paste ("& \\multicolumn{2}{c}{\\emph{winner}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(federal)}} \\\\ \\cmidrule{2-3} \\cmidrule{5-6} \\cmidrule{8-9} \\cmidrule{11-12} \n")
Header3 <- paste ("\\multicolumn{1}{l}{(a) \\textsc{pj}} & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{12}{l}{(b) \\textsc{ucr}} \\\\ \\midrule \n")

addtorow <- list()
addtorow$pos <- list()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 6
addtorow$command <- c (Header1, Header2, Header3, Header4)
print (xtable ( cbind (rows, tabMain)
         , align=c("l","l","c","c","c","c","c","c","c","c","c","c","c")
         , digits=2
         , caption="{Mayoral incumbency effects in the province of Buenos Aires, 1983-2015}"
         , label="T:main")
    , sanitize.text.function=function(x){x}
    , floating=TRUE
    , table.placement="t"
    , caption.placement="top" 
    , latex.environments="center"
    , size="footnotesize"
    , include.colnames=FALSE
    , include.rownames=FALSE
    , hline.after = c ()
    , add.to.row=addtorow )


## (3.1.2) heterogeneous effects (I): by copartisanship status (Table A10) ####
tabCopartisan <- rbind (
  matrix (rd_results_copartisan[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),2], byrow=F, nrow=6)
  , matrix (rd_results_copartisan[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),3], byrow=F, nrow=6)
  , matrix (rd_results_copartisan[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),6], byrow=F, nrow=6)
  , matrix (rd_results_copartisan[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),7], byrow=F, nrow=6) )
(tabCopartisan <- cbind (tabCopartisan[,1:2], rep ("", nrow (tabCopartisan)), tabCopartisan[,3:4], rep ("", nrow (tabCopartisan)), tabCopartisan[,5:6], rep ("", nrow (tabCopartisan)), tabCopartisan[,7:8]))

rows <- rep ( c ("estimate ($\\hat{\\tau}_{\\textsc{rd}}$)"
                 , "95\\% \\textsc{ci}"
                 , "$p$-value"
                 , "bwd."
                 , "$N$"
                 , "control mean"), 4)

Header1 <- paste ("\\toprule & \\multicolumn{2}{c}{} & & \\multicolumn{2}{c}{\\emph{vote share}} & & \\multicolumn{2}{c}{\\emph{seat share}} & & \\multicolumn{2}{c}{\\emph{vote share}} \\\\ \n")
Header2 <- paste ("& \\multicolumn{2}{c}{\\emph{winner}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(federal)}} \\\\ \\cmidrule{2-3} \\cmidrule{5-6} \\cmidrule{8-9} \\cmidrule{11-12} \n")
Header3 <- paste ("\\multicolumn{1}{l}{(a) \\textsc{pj}, copartisan} & \\multicolumn{1}{c}{$t-2$} & \\multicolumn{1}{c}{$t-4$} & & \\multicolumn{1}{c}{$t-2$} & \\multicolumn{1}{c}{$t-4$} & & \\multicolumn{1}{c}{$t-2$} & \\multicolumn{1}{c}{$t-4$} & & \\multicolumn{1}{c}{$t-2$} & \\multicolumn{1}{c}{$t-4$} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{12}{l}{(b) \\textsc{pj}, opposition} \\\\ \\midrule \n")
Header5 <- paste ("[2.5ex] \\multicolumn{12}{l}{(c) \\textsc{ucr}, copartisan} \\\\ \\midrule \n")
Header6 <- paste ("[2.5ex] \\multicolumn{12}{l}{(d) \\textsc{ucr}, opposition} \\\\ \\midrule \n")

addtorow <- list()
addtorow$pos <- list()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 6
addtorow$pos[[5]] <- 12
addtorow$pos[[6]] <- 18
addtorow$command <- c (Header1, Header2, Header3, Header4, Header5, Header6)
print (xtable ( cbind (rows, tabCopartisan)
                , align=c("l","l","c","c","c","c","c","c","c","c","c","c","c")
                , digits=2
                , caption="Robustness checks (\\textsc{v}): Heterogeneous effects by president copartisanship"
                , label="T:robCopartisan")
       , sanitize.text.function=function(x){x}
       , floating=TRUE
       , table.placement="t"
       , caption.placement="top"
       , latex.environments="center"
       , size="footnotesize"
       , include.colnames=FALSE
       , include.rownames=FALSE
       , hline.after = c ()
       , add.to.row=addtorow )


## (3.1.3) heterogeneous effects (II): by net approval (Table A11) ####
tabFavored <- rbind (
  matrix (rd_results_favored[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),2], byrow=F, nrow=6)
  , matrix (rd_results_favored[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),3], byrow=F, nrow=6) )
(tabFavored <- cbind (tabFavored[,1:2], rep ("", nrow (tabFavored)), tabFavored[,3:4], rep ("", nrow (tabFavored)), tabFavored[,5:6], rep ("", nrow (tabFavored)), tabFavored[,7:8]))

rows <- rep ( c ("estimate ($\\hat{\\tau}_{\\textsc{rd}}$)"
                 , "95\\% \\textsc{ci}"
                 , "$p$-value"
                 , "bwd."
                 , "$N$"
                 , "control mean"), 2)

Header1 <- paste ("\\toprule & \\multicolumn{2}{c}{} & & \\multicolumn{2}{c}{\\emph{vote share}} & & \\multicolumn{2}{c}{\\emph{seat share}} & & \\multicolumn{2}{c}{\\emph{vote share}} \\\\ \n")
Header2 <- paste ("& \\multicolumn{2}{c}{\\emph{winner}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(federal)}} \\\\ \\cmidrule{2-3} \\cmidrule{5-6} \\cmidrule{8-9} \\cmidrule{11-12} \n")
Header3 <- paste ("\\multicolumn{1}{l}{(a) High approval} & \\multicolumn{1}{c}{$t-2$} & \\multicolumn{1}{c}{$t-4$} & & \\multicolumn{1}{c}{$t-2$} & \\multicolumn{1}{c}{$t-4$} & & \\multicolumn{1}{c}{$t-2$} & \\multicolumn{1}{c}{$t-4$} & & \\multicolumn{1}{c}{$t-2$} & \\multicolumn{1}{c}{$t-4$} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{12}{l}{(b) Low approval} \\\\ \\midrule \n")

addtorow <- list()
addtorow$pos <- list()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 6
addtorow$command <- c (Header1, Header2, Header3, Header4)
print (xtable ( cbind (rows, tabFavored)
                , align=c("l","l","c","c","c","c","c","c","c","c","c","c","c")
                , digits=2
                , caption="Robustness checks (\\textsc{vi}): Heterogeneous effects by presidential approval (\\textsc{pj} only)"
                , label="T:robApproval")
       , sanitize.text.function=function(x){x}
       , floating=TRUE
       , table.placement="t"
       , caption.placement="top"
       , latex.environments="center"
       , size="footnotesize"
       , include.colnames=FALSE
       , include.rownames=FALSE
       , hline.after = c ()
       , add.to.row=addtorow )





#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
#### // 4 // ADDITIONAL TABLES & FIGURES ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#

### (4.1) Descriptive statistics ####

## (4.1.1) Table w/descriptive statistics (Table A3) ####

# full sample sizes, PJ and UCR
filter (pba, year < 2019) %>% group_by (party, midterm) %>% summarise (N = n())

# we first select the variables to report and use pivot_longer() to create a "long" table
tabDescriptive <- select (filter (pba, !is.na (running) & year < 2019), party, midterm
        , running, winner, voteSh
        , all_of (out_main)
        , copartisan_presi_t2, copartisan_presi_t4, approval_high_t2, approval_high_t4
        ) %>% mutate (
         winner = winner*100
         , approval_high_t2 = ifelse (copartisan_presi_t2==0, NA, approval_high_t2)
         , approval_high_t4 = ifelse (copartisan_presi_t4==0, NA, approval_high_t4)
         ) %>% pivot_longer (
          cols = running:approval_high_t4, names_to = "variable")
tabDescriptive <- tabDescriptive %>% 
  filter (!is.na (value)) %>% 
  group_by (party, midterm, variable) %>% 
  summarise (
    N = n()
    , mean = mean (value)
    , sd = sd (value)
    , min = min (value)
    , max = max (value)
    , whitespace = "") %>% ungroup () %>% mutate (
      variable = factor (variable, levels=c (
        "winner", "running", "voteSh"
        , "winner_t2", "winner_t4", "voteSh_t2", "voteSh_t4"
        , "concejalSh_t2", "concejalSh_t4", "dipuSh_t2", "presiSh_t4"
        , "copartisan_presi_t2", "copartisan_presi_t4", "approval_high_t2", "approval_high_t4")) )

# updating variable names
tabDescriptive <- tabDescriptive %>% 
  arrange (party, midterm, variable) %>% 
  mutate (
    varname = as.character (variable)
    , varname = recode (
      varname,
      
      ## explanatory, running and related variable(s)
      "winner" = "\\emph{winner}$_{t}$ (0/100)"
      , "running" = "\\emph{margin of victory}$_{t}$ (-100:100)"
      , "voteSh" = "\\emph{vote share (municipal)}$_{t}$ (0:100)"
      
      ## outcome variables
      , "winner_t2" = "[1.5ex] \\emph{winner}$_{t+2}$ (0/100)"
      , "winner_t4" = "\\emph{winner}$_{t+4}$ (0/100)"
      , "voteSh_t2" = "\\emph{vote share (municipal)}$_{t+2}$ (0:100)"
      , "voteSh_t4" = "\\emph{vote share (municipal)}$_{t+4}$ (0:100)"
      , "concejalSh_t2" = "\\emph{seat share (municipal)}$_{t+2}$ (0:100)"
      , "concejalSh_t4" = "\\emph{seat share (municipal)}$_{t+4}$ (0:100)"
      , "dipuSh_t2" = "\\emph{vote share (national deputy)}$_{t+2}$ (0:100)"
      , "presiSh_t4" = "\\emph{vote share (president)}$_{t+4}$ (0:100)"
      
      ## intervening variables
      , "copartisan_presi_t2" = "[1.5ex] \\emph{copartisan president}$_{t+2}$ (0/1)"
      , "copartisan_presi_t4" = "\\emph{copartisan president}$_{t+4}$ (0/1)"
      , "approval_high_t2" = "\\emph{popular copartisan president}$_{t+2}$ (0/1)"
      , "approval_high_t4" = "\\emph{popular copartisan president}$_{t+4}$ (0/1)"
      ))

# exporting the table to LaTeX
tabDescriptive <- bind_cols (select (filter (tabDescriptive, party=="pj"), varname, N, mean, sd, min, max, whitespace), select (filter (tabDescriptive, party=="ucr"), N, mean, sd, min, max))
tabDescriptive <- bind_cols (
 as.matrix (tabDescriptive[,1])
 , format (round (as.matrix (tabDescriptive[,2]), 0))
 , format (round (as.matrix (tabDescriptive[,3:6]), digits))
 , as.matrix (tabDescriptive[,7])
 , format (round (as.matrix (tabDescriptive[,8]), 0))
 , format (round (as.matrix (tabDescriptive[,9:12]), digits)) )
tabDescriptive

Header1 <- paste ("\\toprule & \\multicolumn{5}{c}{\\textsc{pj}} & & \\multicolumn{5}{c}{\\textsc{ucr}} \\\\ \n")
Header2 <- paste ("& \\multicolumn{5}{c}{($\\textsc{n} = ", sprintf ("%.0f", nrow (filter (pba, !is.na (running) & party=="pj" & midterm=="concurrent" & year<2019))), "$)} & & \\multicolumn{5}{c}{($\\textsc{n} = ", sprintf ("%.0f", nrow (filter (pba, !is.na (running) & party=="ucr" & midterm=="concurrent" & year<2019))), "$)} \\\\ \\cmidrule{2-6} \\cmidrule{8-12} \n")
Header3 <- paste ("\\multicolumn{1}{l}{(a) Concurrent election years} & \\multicolumn{1}{c}{$N$} & \\multicolumn{1}{c}{mean} & \\multicolumn{1}{c}{sd.} & \\multicolumn{1}{c}{min} & \\multicolumn{1}{c}{max} & & \\multicolumn{1}{c}{$N$} & \\multicolumn{1}{c}{mean} & \\multicolumn{1}{c}{sd.} & \\multicolumn{1}{c}{min} & \\multicolumn{1}{c}{max} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{1}{l}{(b) Midterm election years} & \\multicolumn{5}{c}{($\\textsc{n} = ", sprintf ("%.0f", nrow (filter (pba, !is.na (running) & party=="pj" & midterm=="midterm" & year<2019))), "$)} & & \\multicolumn{5}{c}{($\\textsc{n} = ", sprintf ("%.0f", nrow (filter (pba, !is.na (running) & party=="ucr" & midterm=="midterm" & year<2019))), "$)} \\\\ \\midrule \n")

addtorow <- list()
addtorow$pos <- list()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 15
addtorow$command <- c (Header1, Header2, Header3, Header4)
print (xtable ( as.matrix (tabDescriptive)
         , align=c("l","l","r","r","r","r","r","r","r","r","r","r","r")
         , digits=2
         , caption="{Descriptive statistics}"
         , label="T:descriptive")
    , sanitize.text.function=function(x){x}
    , floating=TRUE
    , table.placement="t"
    , caption.placement="top" 
    , latex.environments="center"
    , size="footnotesize"
    , include.colnames=FALSE
    , include.rownames=FALSE
    , hline.after = c ()
    , add.to.row=addtorow )


## (4.1.2) Correlations between outcomes ####

## (4.1.2.1) correlations

# computing the correlations
out_labels <- c ("winner (t+2)", "winner (t+4)", "vote share (t+2)\n(municipal)", "vote share (t+4)\n(municipal)"
         , "seat share (t+2)\n(municipal)", "seat share (t+4)\n(municipal)", "vote share (t+2)\n(national deputy)", "vote share (t+4)\n(president)")

pba <- pba %>% ungroup () ## we may get an error message otherwise
corr_pj <- round (cor (select (filter (pba, !is.na (running) & year < 2019 & party=="pj"), all_of (out_main)), use="pairwise.complete.obs"), digits)
row.names (corr_pj) <- out_labels; colnames (corr_pj) <- out_labels
(corr_pj_pval <- round (cor_pmat (select (filter (pba, !is.na (running) & year < 2019 & party=="pj"), all_of (out_main))), digits+1)) ## only p-value greater than 0.01 is the correlation between concejalMaj_t2 and presiSh_t4 (p = 0.113)

corr_ucr <- round (cor (select (filter (pba, !is.na (running) & year < 2019 & party=="ucr"), all_of (out_main)), use="pairwise.complete.obs"), digits)
row.names (corr_ucr) <- out_labels; colnames (corr_ucr) <- out_labels
(corr_ucr_pval <- round (cor_pmat (select (filter (pba, !is.na (running) & year < 2019 & party=="ucr"), all_of (out_main))), digits+1)) ## EVERYTHING is significant at the 0.01 level

# drawing the plots: Figure A1a and A1b
(pcorr_pj <- ggcorrplot (corr_pj, method = "square", type = "lower"
              , lab=T, lab_size=size_text*1.5, colors = c (col_density_left, "white", col_density_right) ))
(pcorr_ucr <- ggcorrplot (corr_ucr, method = "square", type = "lower"
              , lab=T, lab_size=size_text*1.5, colors = c (col_density_left, "white", col_density_right) ))


## (4.1.2.2) visualizing the relationship between municipal and federal vote shares
dpj <- dpj %>% ungroup (); ducr <- ducr %>% ungroup ()
range (select (filter (dpj, !is.na (running)), voteSh, fedSh), na.rm=T)
range (select (filter (ducr, !is.na (running)), voteSh, fedSh), na.rm=T) ## between 0 and 80

# Figure A2a
(pcorrsh_pj <- ggplot (filter (dpj, !is.na (running) & year < 2019), aes (x=voteSh, y=fedSh, col=midterm))
 + geom_abline (linetype=2, col=col_concurrent)
 + geom_smooth (method="lm", se=F)
 + geom_point (alpha=alpha_values*3, size=size_bin)
 + scale_color_manual (name="", labels=c ("concurrent years", "midterm years"), values=c (col_line, col_bin))
 + xlim (0, 80) + ylim (0, 80)
 + xlab (expression (vote~share~(municipal)[t]~('%')))
 + ylab (expression (vote~share~(federal)[t]~('%')))
 + theme (legend.position="bottom", legend.title=element_blank (), legend.box.margin=margin (-15,-9,-9,-9)) )

# Figure A2b
(pcorrsh_ucr <- ggplot (filter (ducr, !is.na (running) & year < 2019), aes (x=voteSh, y=fedSh, col=midterm))
 + geom_abline (linetype=2, col=col_concurrent)
 + geom_smooth (method="lm", se=F)
 + geom_point (alpha=alpha_values*3, size=size_bin)
 + scale_color_manual (name="", labels=c ("concurrent years", "midterm years"), values=c (col_line, col_bin))
 + xlim (0, 80) + ylim (0, 80)
 + xlab (expression (vote~share~(municipal)[t]~('%')))
 + ylab (expression (vote~share~(federal)[t]~('%')))
 + theme (legend.position="bottom", legend.title=element_blank (), legend.box.margin=margin (-15,-9,-9,-9)) )

# all estimates are statistically significant at the 0.01 level:
summary (felm (fedSh ~ voteSh*midterm, data=filter (dpj, !is.na (running) & year < 2019)))
summary (felm (fedSh ~ voteSh*midterm, data=filter (ducr, !is.na (running) & year < 2019)))


## (4.1.2.3) exporting the plots
pwid <- 15*1.05
phei <- 15*1.05

ggsave ("figCorrelPJ.png"
    , pcorr_pj, width=pwid, height=phei, units="cm")
ggsave ("figCorrelUCR.png"
    , pcorr_ucr, width=pwid, height=phei, units="cm")

pwid <- 10*1.05
phei <- 10*1.05

ggsave ("figCorrelShPJ.png"
    , pcorrsh_pj, width=pwid, height=phei, units="cm")
ggsave ("figCorrelShUCR.png"
    , pcorrsh_ucr, width=pwid, height=phei, units="cm")


## (4.1.3) Time series (Figure A3) ####

# municipal victories
(pWinner <- ggplot (pba %>% group_by (year, party) %>% summarise (winner=mean (winner, na.rm=T)*100)
          , aes (x=year, y=winner, col=party, linetype=party))
 + geom_vline (xintercept=unique (filter (pba, midterm=="concurrent")$year), col=col_concurrent)
 + geom_hline (yintercept=50, linetype=2, col=col_cutoff)
 + geom_line (size=size_tsline)
 + geom_jitter (data=pba, aes (x=year_party, y=winner*100, col=party), alpha=alpha_values, width=jitter_w, height=jitter_h)
 + xlim (1983, 2019) + ylim (0, 100)
 + scale_x_continuous (breaks=unique (filter (pba, midterm=="concurrent")$year))
 + scale_color_manual (name="", labels=c ("PJ", "UCR"), values=c(col_pj, col_ucr))
 + scale_linetype_discrete (name="", labels=c ("PJ", "UCR"))
 + xlab ("") + ylab ("municipal\nvictories (%)") )

# vote share in municipal election
(pVShare <- ggplot (pba %>% group_by (year, party) %>% summarise (voteSh=mean (voteSh, na.rm=T))
          , aes (x=year, y=voteSh, col=party, linetype=party))
 + geom_vline (xintercept=unique (filter (pba, midterm=="concurrent")$year), col=col_concurrent)
 + geom_hline (yintercept=50, linetype=2, col=col_cutoff)
 + geom_line (size=size_tsline)
 + geom_jitter (data=pba, aes (x=year_party, y=voteSh, col=party), alpha=alpha_values, width=jitter_w, height=0)
 + xlim (1983, 2019) + ylim (0, 100)
 + scale_x_continuous (breaks=unique (filter (pba, midterm=="concurrent")$year))
 + scale_color_manual (name="", labels=c ("PJ", "UCR"), values=c(col_pj, col_ucr))
 + scale_linetype_discrete (name="", labels=c ("PJ", "UCR"))
 + xlab ("") + ylab ("vote in\nmunicipal election (%)") )

# seat share in municipal election
(pSShare <- ggplot (pba %>% group_by (year, party) %>% summarise (concejalSh=mean (concejalSh, na.rm=T))
          , aes (x=year, y=concejalSh, col=party, linetype=party))
 + geom_vline (xintercept=unique (filter (pba, midterm=="concurrent")$year), col=col_concurrent)
 + geom_hline (yintercept=50, linetype=2, col=col_cutoff)
 + geom_line (size=size_tsline)
 + geom_jitter (data=pba, aes (x=year_party, y=concejalSh, col=party), alpha=alpha_values, width=jitter_w, height=0)
 + xlim (1983, 2019) + ylim (0, 100)
 + scale_x_continuous (breaks=unique (filter (pba, midterm=="concurrent")$year))
 + scale_color_manual (name="", labels=c ("PJ", "UCR"), values=c(col_pj, col_ucr))
 + scale_linetype_discrete (name="", labels=c ("PJ", "UCR"))
 + xlab ("") + ylab ("seats captured in\nmunicipal election (%)") )

# vote share in federal election
(pFedShare <- ggplot (pba %>% group_by (year, party) %>% summarise (fedSh=mean (fedSh, na.rm=T))
          , aes (x=year, y=fedSh, col=party, linetype=party))
 + geom_hline (yintercept=50, linetype=2, col=col_cutoff)
 + geom_vline (xintercept=unique (filter (pba, year_presid==1)$year), col=col_concurrent)
 + geom_line (size=size_tsline)
 + geom_jitter (data=pba, aes (x=year_party, y=fedSh, col=party), alpha=alpha_values, width=jitter_w, height=0)
 + xlim (1983, 2019) + ylim (0, 80)
 + scale_x_continuous (breaks=unique (filter (pba, year_presid==1)$year))
 + scale_color_manual (name="", labels=c ("PJ", "UCR"), values=c(col_pj, col_ucr))
 + scale_linetype_discrete (name="", labels=c ("PJ", "UCR"))
 + xlab ("") + ylab ("vote in\nfederal election (%)") )


## (4.1.4) Distribution of running variable by seccion (Figure A5) ####

# margins of victory
range (filter (pba, midterm=="concurrent" & year <= 2015)$running, na.rm=TRUE)
(pRunning <- ggplot (pba %>% filter (midterm=="concurrent" & year <= 2015) %>% group_by (seccion, party) %>% summarise (running=mean (running, na.rm=T)) %>% mutate (seccion=as.numeric (seccion))
          , aes (x=seccion, y=running, col=party, linetype=party))
 + annotate ("rect", xmin = 0.5, xmax = 8.5, ymin = -max (as.numeric (tabMain[c (4, 10),]), na.rm=TRUE), ymax = max (as.numeric (tabMain[c (4, 10),]), na.rm=TRUE), alpha = alpha_values*2)
 + annotate ("rect", xmin = 0.5, xmax = 8.5, ymin = -median (as.numeric (tabMain[c (4, 10),]), na.rm=TRUE), ymax = median (as.numeric (tabMain[c (4, 10),]), na.rm=TRUE), alpha = alpha_values*2)
 + geom_hline (yintercept=0, linetype=2, col=col_cutoff)
 + geom_line (size=size_tsline)
 + geom_jitter (data=pba %>% filter (midterm=="concurrent"), aes (x=seccion_party, y=running, col=party, shape=conurba), alpha=alpha_values*2, width=jitter_w/2, height=jitter_h/2)
 + scale_y_continuous (breaks=seq (-60, 60, by=15))
 + scale_x_continuous (limits=c (0.5, 8.5), breaks=1:8, labels=str_c (levels (pba$seccion), c ("*", "", "*", rep ("", 4), "*")))
 + scale_color_manual (name="", labels=c ("PJ", "UCR"), values=c(col_pj, col_ucr))
 + scale_shape_discrete (name="")
 + scale_linetype_discrete (name="", labels=c ("PJ", "UCR"))
 + xlab ("electoral 'sección'") + ylab ("margin in\nmayoral election (%)")
  )


## exporting the plots
pwid <- 20*1.325
phei <- 20*0.375

ggsave ("figTimeSeriesWinner.png"
    , pWinner, width=pwid, height=phei, units="cm")
ggsave ("figTimeSeriesVShare.png"
    , pVShare, width=pwid, height=phei, units="cm")
ggsave ("figTimeSeriesSeats.png"
    , pSShare, width=pwid, height=phei, units="cm")
ggsave ("figTimeSeriesFShare.png"
    , pFedShare, width=pwid, height=phei, units="cm")
ggsave ("figRunning.png"
    , pRunning, width=pwid, height=phei*1.5, units="cm")


## (4.1.5) Maps ####

# downloading the shapefiles
setwd (home)
map <- st_read ("maps/departamento.shp")
sort (unique (map$SAG)) ## we're only interested in "ARBA - Gerencia de Servicios Catastrales
map <- filter (map, SAG=="ARBA - Gerencia de Servicios Catastrales") %>% 
  mutate (
 municipio = toupper (FNA) )

# updating crs
st_crs (map)
map <- st_transform (map, 5340) ## POSGAR 2007 projection

# municipalities that compose the Conurbano
conu <- read_csv ("cordones.csv") %>% 
  mutate (
    municipio = toupper (municipio) )
summary (sort (unique (conu$municipio)) %in% sort (unique (pba$municipio))) ## only La Matanza (Oeste), which is not a municipality but part of one. That's OK:
sort (unique (conu$municipio))[(sort (unique (conu$municipio)) %in% sort (unique (pba$municipio)))==F]
setwd (str_c (home, "figures/"))

# updating municipality names
summary (sort (unique (map$municipio)) %in% sort (unique (pba$municipio))) ## 40 discrepancies
sort (unique (map$municipio))[(sort (unique (map$municipio)) %in% sort (unique (pba$municipio)))==F]
map <- map %>% 
  mutate (
    municipio = str_replace_all (municipio, c (
      "Á" = "A"
      , "É" = "E"
      , "Í" = "I"
      , "Ó" = "O"
      , "Ú" = "U"
      , "Ñ" = "N"
      , "25" = "VEINTICINCO"
      , "9" = "NUEVE"
      , "DE MARINA LEONARDO " = ""
      , "JOSE M. " = ""
      , "SAN CARLOS DE " = ""
      , "LA COSTA" = "PARTIDO DE LA COSTA"
      , "GENERAL JUAN " = "GENERAL "
      )))
summary (sort (unique (map$municipio)) %in% sort (unique (pba$municipio))) ## no more discrepancies. Good

# joining
map <- left_join (map, conu, by=c ("municipio" = "municipio")) %>% 
  mutate (
 conurbano06 = ifelse (is.na (conurbano06), 0, conurbano06) )
map <- left_join (select (map, municipio, conurbano06, geometry)
         , pba %>% 
           group_by (municipio, party) %>% 
           summarise (
             winner = mean (winner, na.rm=T)*100
             , voteSh = mean (voteSh, na.rm=T)
             , concejalMaj = mean (concejalMaj, na.rm=T)*100
             , fedSh = mean (fedSh, na.rm=T) )
         , by=c ("municipio" = "municipio")
         )
summary (map)


## Conurbano only

# PJ, municipal vote share (Figure 1a)
(mapConuVShare_pj <- ggplot (filter (map, party=="pj" & conurbano06==1))
 + geom_sf (aes (fill=voteSh), color=col_concurrent, size=alpha_values)
 + scale_fill_viridis (name="", option=map_palette, direction=-1, limits=c (10, 60))
 + theme_map () )

# UCR, municipal vote share (Figure 1b)
(mapConuVShare_ucr <- ggplot (filter (map, party=="ucr" & conurbano06==1))
 + geom_sf (aes (fill=voteSh), color=col_concurrent, size=alpha_values)
 + scale_fill_viridis (name="", option=map_palette, direction=-1, limits=c (10, 60))
 + theme_map () )

# PJ, federal vote share (Figure A4a)
(mapConuFedShare_pj <- ggplot (filter (map, party=="pj" & conurbano06==1))
 + geom_sf (aes (fill=fedSh), color=col_concurrent, size=alpha_values)
 + scale_fill_viridis (name="", option=map_palette, direction=-1, limits=c (10, 60))
 + theme_map () )

# UCR, federal vote share (Figure A4b)
(mapConuFedShare_ucr <- ggplot (filter (map, party=="ucr" & conurbano06==1))
 + geom_sf (aes (fill=fedSh), color=col_concurrent, size=alpha_values)
 + scale_fill_viridis (name="", option=map_palette, direction=-1, limits=c (10, 60))
 + theme_map () )


## non-Conurbano only

# PJ, municipal vote share (Figure 1c)
(mapRestVShare_pj <- ggplot (filter (map, party=="pj" & conurbano06==0))
 + geom_sf (aes (fill=voteSh), color=col_concurrent, size=alpha_values)
 + scale_fill_viridis (name="", option=map_palette, direction=-1, limits=c (10, 60))
 + theme_map () )

# UCR, municipal vote share (Figure 1d)
(mapRestVShare_ucr <- ggplot (filter (map, party=="ucr" & conurbano06==0))
 + geom_sf (aes (fill=voteSh), color=col_concurrent, size=alpha_values)
 + scale_fill_viridis (name="", option=map_palette, direction=-1, limits=c (10, 60))
 + theme_map () )

# PJ, federal vote share (Figure A4c)
(mapRestFedShare_pj <- ggplot (filter (map, party=="pj" & conurbano06==0))
 + geom_sf (aes (fill=fedSh), color=col_concurrent, size=alpha_values)
 + scale_fill_viridis (name="", option=map_palette, direction=-1, limits=c (10, 60))
 + theme_map () )

# UCR, federal vote share (Figure A4d)
(mapRestFedShare_ucr <- ggplot (filter (map, party=="ucr" & conurbano06==0))
 + geom_sf (aes (fill=fedSh), color=col_concurrent, size=alpha_values)
 + scale_fill_viridis (name="", option=map_palette, direction=-1, limits=c (10, 60))
 + theme_map () )

# exporting the plots
pwid <- 14*1.05
phei <- 14*1.15

ggsave ("figMapConuPJmuni.png"
    , mapConuVShare_pj, width=pwid, height=phei, units="cm")
ggsave ("figMapConuUCRmuni.png"
    , mapConuVShare_ucr, width=pwid, height=phei, units="cm")
ggsave ("figMapConuPJfed.png"
    , mapConuFedShare_pj, width=pwid, height=phei, units="cm")
ggsave ("figMapConuUCRfed.png"
    , mapConuFedShare_ucr, width=pwid, height=phei, units="cm")
ggsave ("figMapRestPJmuni.png"
    , mapRestVShare_pj, width=pwid, height=phei, units="cm")
ggsave ("figMapRestUCRmuni.png"
    , mapRestVShare_ucr, width=pwid, height=phei, units="cm")
ggsave ("figMapRestPJfed.png"
    , mapRestFedShare_pj, width=pwid, height=phei, units="cm")
ggsave ("figMapRestUCRfed.png"
    , mapRestFedShare_ucr, width=pwid, height=phei, units="cm")


## (4.1.6) Table A1: Copartisanship and presidential approval ####
(tabPresi <- left_join (copartisan
                        , ead %>% filter (election==1) %>% select (year, approval, approval_high)
                        , by="year") %>% 
   filter (year %in% 1985:2019) %>% 
   mutate (
     eltype = ifelse (
       year %in% seq (1985, 2025, by=4), "midterm", "concurrent")
     , year = as.character (year)
     , presi_pj = ifelse (party_presi=="pj", "\\textsc{pj}", "")
     , presi_ucr = ifelse (party_presi=="ucr", "\\textsc{ucr}", "")
     , gover_pj = ifelse (party_gover=="pj", "\\textsc{pj}", "")
     , gover_ucr = ifelse (party_gover=="ucr", "\\textsc{ucr}", "")
     , net_approval = ifelse (
       is.na (approval), "", ifelse (
         approval_high==0
         , str_c ("{\\color{red} ", sprintf ("%.1f", approval), "}")
         , sprintf ("%.1f", approval)))) %>% 
   select (year, eltype:gover_ucr, net_approval))

Header1 <- paste ("\\toprule \\multicolumn{2}{c}{outcome measured in} & & \\multicolumn{2}{c}{president's party} & & \\multicolumn{2}{c}{governor's party} & & \\multicolumn{1}{c}{net approval \\%} \\\\ \\cmidrule{1-2} \\cmidrule{4-5} \\cmidrule{7-8} \\cmidrule{10-10} \n")

addtorow <- list()
addtorow$pos <- list()
addtorow$pos[[1]] <- 0
addtorow$command <- c (Header1)
print (xtable ( cbind (as.matrix (tabPresi)[,1:2], rep ("", nrow (tabPresi)), as.matrix (tabPresi)[,3:4], rep ("", nrow (tabPresi)), as.matrix (tabPresi)[,5:6], rep ("", nrow (tabPresi)), as.matrix (tabPresi)[,7])
                , align=c("l","l","c","c","c","c","c","c","c","c","r")
                , digits=2
                , caption="Copartisanship and presidential approval, 1985-2019"
                , label="T:presidents")
       , sanitize.text.function=function(x){x}
       , floating=TRUE
       , table.placement="t"
       , caption.placement="top" 
       , latex.environments="center"
       , size="footnotesize"
       , include.colnames=FALSE
       , include.rownames=FALSE
       , hline.after = c ()
       , add.to.row=addtorow )



### (4.2) Balance checks ####

## (4.2.1) RD density plots ####

## PJ, concurrent elections (Figure A8a)
density_pj_conc <- rddensity (filter (dpj, !is.na(running) & midterm=="concurrent" & year<=2015)$running) ## we exclude 2019 b/c we have no data for future elections
pDensityPJ_conc <- rdplotdensity (rdd=density_pj_conc
                 , X=filter (dpj, !is.na(running) & midterm=="concurrent" & year<=2015)$running
                 , type="both", plotRange=c(-50, 50)
                 , lcol=c (col_density_left, col_density_right)
                 , CIcol=c (col_density_left, col_density_right), CIshade=shade_density_ci
                 , histFillCol=col_density_hist, histFillShade=shade_density_hist)
(pDensityPJ_conc <- pDensityPJ_conc$Estplot
 + annotate ("label", label=paste0 ("p-value: ", format (round (density_pj_conc$test$p_jk, digits), nsmall=digits)), x=35, y=0.04, size=size_text)
 + theme (text=element_text(size=16)) )

# UCR, concurrent elections (Figure A8b)
density_ucr_conc <- rddensity (filter (ducr, !is.na(running) & midterm=="concurrent" & year<=2015)$running)
pDensityUCR_conc <- rdplotdensity (rdd=density_ucr_conc
                  , X=filter (ducr, !is.na(running) & midterm=="concurrent" & year<=2015)$running
                  , type="both", plotRange=c(-50, 50)
                  , lcol=c (col_density_left, col_density_right)
                  , CIcol=c (col_density_left, col_density_right), CIshade=shade_density_ci
                  , histFillCol=col_density_hist, histFillShade=shade_density_hist)
(pDensityUCR_conc <- pDensityUCR_conc$Estplot
 + annotate ("label", label=paste0 ("p-value: ", format (round (density_ucr_conc$test$p_jk, digits), nsmall=digits)), x=35, y=0.045, size=size_text)
 + theme (text=element_text(size=16)) )

# PJ, midterm elections (Figure A8c)
density_pj_midt <- rddensity (filter (dpj, !is.na(running) & midterm=="midterm" & year<=2017)$running)
pDensityPJ_midt <- rdplotdensity (rdd=density_pj_midt
                 , X=filter (dpj, !is.na(running) & midterm=="midterm" & year<=2017)$running
                 , type="both", plotRange=c(-50, 50)
                 , lcol=c (col_density_left, col_density_right)
                 , CIcol=c (col_density_left, col_density_right), CIshade=shade_density_ci
                 , histFillCol=col_density_hist, histFillShade=shade_density_hist)
(pDensityPJ_midt <- pDensityPJ_midt$Estplot
 + annotate ("label", label=paste0 ("p-value: ", format (round (density_pj_midt$test$p_jk, digits), nsmall=digits)), x=35, y=0.03, size=size_text)
 + theme (text=element_text(size=16)) )

# UCR, midterm elections (Figure A8d)
density_ucr_midt <- rddensity (X = filter (ducr, !is.na(running) & midterm=="midterm" & year<=2017)$running)
pDensityUCR_midt <- rdplotdensity (rdd=density_ucr_midt
                  , X=filter (ducr, !is.na(running) & midterm=="midterm" & year<=2017)$running
                  , type="both", plotRange=c(-50, 50)
                  , lcol=c (col_density_left, col_density_right)
                  , CIcol=c (col_density_left, col_density_right), CIshade=shade_density_ci
                  , histFillCol=col_density_hist, histFillShade=shade_density_hist)
(pDensityUCR_midt <- pDensityUCR_midt$Estplot
 + annotate ("label", label=paste0 ("p-value: ", format (round (density_ucr_midt$test$p_jk, digits), nsmall=digits)), x=35, y=0.03, size=size_text)
 + theme (text=element_text(size=16)) )

# getting the p-values
round (density_pj_conc$test$p_jk, 2) ## p-value: 0.93
round (density_ucr_conc$test$p_jk, 2) ## p-value: 0.56
round (density_pj_midt$test$p_jk, 2) ## p-value: 0.80
round (density_ucr_midt$test$p_jk, 2) ## p-value: 0.65

# exporting the plots
pwid <- 12.5*1.05
phei <- 12.5*0.90

ggsave ("figDensityPJconc.png"
    , pDensityPJ_conc, width=pwid, height=phei, units="cm")
ggsave ("figDensityUCRconc.png"
    , pDensityUCR_conc, width=pwid, height=phei, units="cm")
ggsave ("figDensityPJmidt.png"
    , pDensityPJ_midt, width=pwid, height=phei, units="cm")
ggsave ("figDensityUCRmidt.png"
    , pDensityUCR_midt, width=pwid, height=phei, units="cm")


## (4.2.2) lagged outcomes ####

## doing the analysis
rd_results_copartisan_lag <- NULL; rd_results_favored_lag <- NULL
for (i in 1:length (out_lag)){
  
  ## PJ
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_lag[i])])
  
  all <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_lag2", out_lag)) { ## outcome measured at t-2
    hi <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_lag2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_lag2==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_lag4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_lag4==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  }
  
  pj_copartisan_lag <- rd_hte (all, hi, lo, out_lag[i], "PJ", "Copartisan", "Opposition"
                               , mean_y_all = mean (filter (dpj, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                               , mean_y_hi = ifelse (
                                 i %in% grep ("_lag2", out_lag)
                                 , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_lag2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                 , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_lag4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                               , mean_y_lo = ifelse (
                                 i %in% grep ("_lag2", out_lag)
                                 , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_lag2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                 , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_lag4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## UCR
  ducr$outcome_tmp <- unlist (ducr[,which (colnames (ducr) == out_lag[i])])
  
  all <- with (filter (ducr, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_lag2", out_lag)) { ## outcome measured at t+2
    hi <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_lag2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_lag2==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_lag4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_lag4==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  }
  
  ucr_copartisan_lag <- rd_hte (all, hi, lo, out_lag[i], "UCR", "Copartisan", "Opposition"
                                , mean_y_all = mean (filter (ducr, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                , mean_y_hi = ifelse (
                                  i %in% grep ("_lag2", out_lag)
                                  , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_lag2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                  , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_lag4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                                , mean_y_lo = ifelse (
                                  i %in% grep ("_lag2", out_lag)
                                  , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_lag2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                  , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_lag4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_copartisan_lag <- rbind (rd_results_copartisan_lag, cbind (matrix (unlist (
    list (pj_copartisan_lag, ucr_copartisan_lag)), nrow=8)))
}
rd_results_copartisan_lag


## favored copartisan president: PJ only
for (i in 1:length (out_lag)){
  
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_lag[i])])
  
  all <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_lag2", out_lag)) { ## outcome measured at t+2
    hi <- with (filter (dpj, midterm=="concurrent" & favored_presi_lag2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & disfavored_presi_lag2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (dpj, midterm=="concurrent" & favored_presi_lag4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & disfavored_presi_lag4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  }
  
  pj_favored_lag <- rd_hte (all, hi, lo, out_lag[i], "PJ", "Favored", "Disfavored"
                            , mean_y_all = mean (filter (dpj, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                            , mean_y_hi = ifelse (
                              i %in% grep ("_lag2", out_lag)
                              , mean (filter (dpj, midterm=="concurrent" & favored_presi_lag2==1 & !is.na (favored_presi_lag2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                              , mean (filter (dpj, midterm=="concurrent" & favored_presi_lag4==1 & !is.na (favored_presi_lag4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                            , mean_y_lo = ifelse (
                              i %in% grep ("_lag2", out_lag)
                              , mean (filter (dpj, midterm=="concurrent" & disfavored_presi_lag2==0 & !is.na (disfavored_presi_lag2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                              , mean (filter (dpj, midterm=="concurrent" & disfavored_presi_lag4==0 & !is.na (disfavored_presi_lag4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_favored_lag <- rbind (rd_results_favored_lag, cbind (matrix (unlist (
    list (pj_favored_lag)), nrow=8)))
}
rd_results_favored_lag


## (4.2.2.1) full sample (Table A4) ####
tabFullLag <- rbind (
 matrix (rd_results_copartisan_lag[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),1]
         , byrow=F, nrow=6)
 , matrix (rd_results_copartisan_lag[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),5]
           , byrow=F, nrow=6))
(tabFullLag <- cbind (tabFullLag[,1:2], rep ("", nrow (tabFullLag)), tabFullLag[,3:4], rep ("", nrow (tabFullLag)), tabFullLag[,5:6], rep ("", nrow (tabFullLag)), tabFullLag[,7:8]))

rows <- rep ( c ("estimate ($\\hat{\\tau}_{\\textsc{rd}}$)"
         , "95\\% \\textsc{ci}"
         , "$p$-value"
         , "bwd."
         , "$N$"
         , "control mean"), 2)

Header1 <- paste ("\\toprule & \\multicolumn{2}{c}{} & & \\multicolumn{2}{c}{\\emph{vote share}} & & \\multicolumn{2}{c}{\\emph{seat share}} & & \\multicolumn{2}{c}{\\emph{vote share}} \\\\ \n")
Header2 <- paste ("& \\multicolumn{2}{c}{\\emph{winner}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(federal)}} \\\\ \\cmidrule{2-3} \\cmidrule{5-6} \\cmidrule{8-9} \\cmidrule{11-12} \n")
Header3 <- paste ("\\multicolumn{1}{l}{(a) \\textsc{pj}} & \\multicolumn{1}{c}{$t-2$} & \\multicolumn{1}{c}{$t-4$} & & \\multicolumn{1}{c}{$t-2$} & \\multicolumn{1}{c}{$t-4$} & & \\multicolumn{1}{c}{$t-2$} & \\multicolumn{1}{c}{$t-4$} & & \\multicolumn{1}{c}{$t-2$} & \\multicolumn{1}{c}{$t-4$} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{12}{l}{(b) \\textsc{ucr}} \\\\ \\midrule \n")

addtorow <- list()
addtorow$pos <- list()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 6
addtorow$command <- c (Header1, Header2, Header3, Header4)
print (xtable ( cbind (rows, tabFullLag)
        , align=c("l","l","c","c","c","c","c","c","c","c","c","c","c")
        , digits=2
        , caption="{Balance checks: Mayoral incumbency effects on \\emph{lagged} outcomes}"
        , label="T:balance")
    , sanitize.text.function=function(x){x}
    , floating=TRUE
    , table.placement="t"
    , caption.placement="top" 
    , latex.environments="center"
    , size="footnotesize"
    , include.colnames=FALSE
    , include.rownames=FALSE
    , hline.after = c ()
    , add.to.row=addtorow )





#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#
#### // 5 // ROBUSTNESS CHECKS ####
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%#

### (5.1) Sensitivity to bandwidth choice ####

## alternative bandwidth values + objects for storage
(bwds <- seq (4, 40, by=2))
pj_bwd_sensitivity <- NULL; ucr_bwd_sensitivity <- NULL

## doing the loops
for (o in 1:length (out_main)){
 
 ### selecting the outcomes of interest
 dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == outcomes[o])])
 ducr$outcome_tmp <- unlist (ducr[,which (colnames (ducr) == outcomes[o])])
 
 
 ### for each outcome, we need the optimal bandwidth, as well as half and double the optimal
 rd_pj_tmp_main <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
 rd_ucr_tmp_main <- with (filter (ducr, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
 bwds_pj_tmp <- c (bwds, rd_pj_tmp_main$bws[1,1], rd_pj_tmp_main$bws[1,1]*1/2, rd_pj_tmp_main$bws[1,1]*2)
 bwds_ucr_tmp <- c (bwds, rd_ucr_tmp_main$bws[1,1], rd_ucr_tmp_main$bws[1,1]*1/2, rd_ucr_tmp_main$bws[1,1]*2)
 
 
 ### looping over all bandwidth choices
 for (b in 1:(length (bwds)+3)){
  rd_pj_tmp <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio, h=bwds_pj_tmp[b], b=bwds_pj_tmp[b]))
  pj_bwd_sensitivity <- rbind (
   pj_bwd_sensitivity, c (o, rd_pj_tmp$bws[1,1], rd_pj_tmp$Estimate[1], rd_pj_tmp$ci[3,], sum (rd_pj_tmp$N_h)) )
  
  rd_ucr_tmp <- with (filter (ducr, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio, h=bwds_ucr_tmp[b], b=bwds_ucr_tmp[b]))
  ucr_bwd_sensitivity <- rbind (
   ucr_bwd_sensitivity, c (o, rd_ucr_tmp$bws[1,1], rd_ucr_tmp$Estimate[1], rd_ucr_tmp$ci[3,], sum (rd_ucr_tmp$N_h))
   )}}


## turning into a data frame for plotting
pj_bwd_sensitivity <- as.data.frame (pj_bwd_sensitivity)
colnames (pj_bwd_sensitivity) <- c ("outcome", "bwd", "estim", "low", "high", "N")
pj_bwd_sensitivity$party <- "pj"
ucr_bwd_sensitivity <- as.data.frame (ucr_bwd_sensitivity)
colnames (ucr_bwd_sensitivity) <- c ("outcome", "bwd", "estim", "low", "high", "N")
ucr_bwd_sensitivity$party <- "ucr"

plot_sensitivity <- bind_rows (pj_bwd_sensitivity, ucr_bwd_sensitivity) %>% mutate (outcome = as.character (outcome))
plot_sensitivity <- left_join (as.data.frame (cbind (1:length (out_main), out_main)), plot_sensitivity, by=c ("V1" = "outcome")) %>% 
 separate (out_main, c ("out_var", "out_time"), remove=F, sep="_") %>% mutate (
  V1 = NULL
  , party = factor (party)
  , out_main = factor (out_main)
  
  # out_var: type of outcome variable (victory, vote share, etc)
  , out_var = recode (out_var
            , winner = "victory (0/100)"
            , voteSh = "vote share\n(municipal) (0:100)"
            , concejalSh = "seat share\n(municipal) (0:100)"
            , concejalMaj = "majority in\ncouncil (0/100)"
            , presiSh = "vote share\n(federal) (0:100)"
            , dipuSh = "vote share\n(federal) (0:100)" )
  , out_var = factor (out_var, levels = c (
   "victory (0/100)", "vote share\n(municipal) (0:100)", "seat share\n(municipal) (0:100)", "majority in\ncouncil (0/100)", "vote share\n(federal) (0:100)"))
  
  # out_var: time at which outcome variable is measured (t+2, t+4, etc)
  , out_time = recode (out_time, t2 = "t + 2", t4 = "t + 4")
  , out_time = factor (out_time)
  
  , bwd_type = ifelse (bwd %in% bwds, "manual: 4 to 40", "1/2x; 1x; or 2x CCT-optimal bandwidth")
  , bwd_type = factor (bwd_type) )
summary (plot_sensitivity)
range (select (plot_sensitivity, low, high)) ## -40.11 to 80.46


## Figure A10
(bwdp11_main_pj <- ggplot (filter (plot_sensitivity, party=="pj"), aes (x=bwd, y=estim, color=bwd_type))
 + geom_hline (yintercept=0, col=col_cutoff)
 + geom_linerange (aes (ymin=low, ymax=high), size=size_bin/2)
 + geom_point (size=size_bin, alpha=alpha_bin)
 + scale_color_manual (values=c (col_line, col_bin))
 + facet_grid (out_time ~ out_var)
 + xlim (0, 40)
 + scale_y_continuous (breaks=seq (-25, 100, by=25), limits=c (-40.15, 80.5))
 + xlab ("bandwidth") + ylab ("RD point estimate + 95% CI")
 + theme (legend.position="bottom", legend.title=element_blank (), legend.box.margin=margin (-15,-9,-9,-9))
 + theme (axis.title.x = element_text (margin=margin (t=10, r=0, b=0, l=0))) )

(bwdp21_main_ucr <- ggplot (filter (plot_sensitivity, party=="ucr"), aes (x=bwd, y=estim, color=bwd_type))
 + geom_hline (yintercept=0, col=col_cutoff)
 + geom_linerange (aes (ymin=low, ymax=high), size=size_bin/2)
 + geom_point (size=size_bin, alpha=alpha_bin)
 + scale_color_manual (values=c (col_line, col_bin))
 + facet_grid (out_time ~ out_var)
 + xlim (0, 40)
 + scale_y_continuous (breaks=seq (-25, 100, by=25), limits=c (-40.15, 80.5))
 + xlab ("bandwidth") + ylab ("RD point estimate + 95% CI")
 + theme (legend.position="bottom", legend.title=element_blank (), legend.box.margin=margin (-15,-9,-9,-9))
 + theme (axis.title.x = element_text (margin=margin (t=10, r=0, b=0, l=0))) )


## exporting the plots
pwid <- 17*1.5
phei <- 17*0.875

ggsave ("figBwdPlotPJ.png"
    , bwdp11_main_pj, width=pwid, height=phei, units="cm")
ggsave ("figBwdPlotUCR.png"
    , bwdp21_main_ucr, width=pwid, height=phei, units="cm")



### (5.2) CER-optimal bandwidths

## (5.2.1) running the models ####
rd_results_copartisan_cer <- NULL; rd_results_favored_cer <- NULL
for (i in 1:length (out_main)){
  
  ## PJ
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_main[i])])
  
  all <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_t2", out_main)) { ## outcome measured at t+2
    hi <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==0), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==0), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio))
  }
  
  pj_copartisan_cer <- rd_hte (all, hi, lo, out_main[i], "PJ", "Copartisan", "Opposition"
                               , mean_y_all = mean (filter (dpj, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                               , mean_y_hi = ifelse (
                                 i %in% grep ("_t2", out_main)
                                 , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                 , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                               , mean_y_lo = ifelse (
                                 i %in% grep ("_t2", out_main)
                                 , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                 , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## UCR
  ducr$outcome_tmp <- unlist (ducr[,which (colnames (ducr) == out_main[i])])
  
  all <- with (filter (ducr, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_t2", out_main)) { ## outcome measured at t+2
    hi <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio));
    lo <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==0), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio));
    lo <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==0), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio))
  }
  
  ucr_copartisan_cer <- rd_hte (all, hi, lo, out_main[i], "UCR", "Copartisan", "Opposition"
                                , mean_y_all = mean (filter (ducr, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                , mean_y_hi = ifelse (
                                  i %in% grep ("_t2", out_main)
                                  , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                  , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                                , mean_y_lo = ifelse (
                                  i %in% grep ("_t2", out_main)
                                  , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                  , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_copartisan_cer <- rbind (rd_results_copartisan_cer, cbind (matrix (unlist (
    list (pj_copartisan_cer, ucr_copartisan_cer)), nrow=8)))
}
rd_results_copartisan_cer


## favored copartisan president: PJ only
for (i in 1:length (out_main)){
  
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_main[i])])
  
  all <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_t2", out_main)) { ## outcome measured at t+2
    hi <- with (filter (dpj, midterm=="concurrent" & favored_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & disfavored_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (dpj, midterm=="concurrent" & favored_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & disfavored_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="cerrd", covs=NULL, cluster=municipio))
  }
  
  pj_favored_cer <- rd_hte (all, hi, lo, out_main[i], "PJ", "Favored", "Disfavored"
                            , mean_y_all = mean (filter (dpj, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                            , mean_y_hi = ifelse (
                              i %in% grep ("_t2", out_main)
                              , mean (filter (dpj, midterm=="concurrent" & favored_presi_t2==1 & !is.na (favored_presi_t2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                              , mean (filter (dpj, midterm=="concurrent" & favored_presi_t4==1 & !is.na (favored_presi_t4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                            , mean_y_lo = ifelse (
                              i %in% grep ("_t2", out_main)
                              , mean (filter (dpj, midterm=="concurrent" & disfavored_presi_t2==0 & !is.na (disfavored_presi_t2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                              , mean (filter (dpj, midterm=="concurrent" & disfavored_presi_t4==0 & !is.na (disfavored_presi_t4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_favored_cer <- rbind (rd_results_favored_cer, cbind (matrix (unlist (
    list (pj_favored_cer)), nrow=8)))
}
rd_results_favored_cer


## (5.2.2) Table A6: full sample ####
tabFullSE <- rbind (
 matrix (rd_results_copartisan_cer[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),1]
         , byrow=F, nrow=6)
 , matrix (rd_results_copartisan_cer[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),5]
           , byrow=F, nrow=6))
(tabFullSE <- cbind (tabFullSE[,1:2], rep ("", nrow (tabFullSE)), tabFullSE[,3:4], rep ("", nrow (tabFullSE)), tabFullSE[,5:6], rep ("", nrow (tabFullSE)), tabFullSE[,7:8]))

rows <- rep ( c ("estimate ($\\hat{\\tau}_{\\textsc{rd}}$)"
         , "95\\% \\textsc{ci}"
         , "$p$-value"
         , "bwd."
         , "$N$"
         , "control mean"), 2)

Header1 <- paste ("\\toprule & \\multicolumn{2}{c}{} & & \\multicolumn{2}{c}{\\emph{vote share}} & & \\multicolumn{2}{c}{\\emph{seat share}} & & \\multicolumn{2}{c}{\\emph{vote share}} \\\\ \n")
Header2 <- paste ("& \\multicolumn{2}{c}{\\emph{winner}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(federal)}} \\\\ \\cmidrule{2-3} \\cmidrule{5-6} \\cmidrule{8-9} \\cmidrule{11-12} \n")
Header3 <- paste ("\\multicolumn{1}{l}{(a) \\textsc{pj}} & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{12}{l}{(b) \\textsc{ucr}} \\\\ \\midrule \n")

addtorow <- list()
addtorow$pos <- list()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 6
addtorow$command <- c (Header1, Header2, Header3, Header4)
print (xtable ( cbind (rows, tabFullSE)
        , align=c("l","l","c","c","c","c","c","c","c","c","c","c","c")
        , digits=2
        , caption="{Robustness checks (\\textsc{i}): \\textsc{cer}-optimal bandwidths}"
        , label="T:robCER")
    , sanitize.text.function=function(x){x}
    , floating=TRUE
    , table.placement="t"
    , caption.placement="top" 
    , latex.environments="center"
    , size="footnotesize"
    , include.colnames=FALSE
    , include.rownames=FALSE
    , hline.after = c ()
    , add.to.row=addtorow )



### (5.3) Clustering SEs by year ####

## (5.3.1) running the models ####
rd_results_copartisan_clustyear <- NULL; rd_results_favored_clustyear <- NULL
for (i in 1:length (out_main)){
  
  ## PJ
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_main[i])])
  
  all <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=factor(year)))
  if (i %in% grep ("_t2", out_main)) { ## outcome measured at t+2
    hi <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=factor(year)));
    lo <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=factor(year)))
  } else {
    hi <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=factor(year)));
    lo <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==0 & !is.na (running)), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL)) ## in this specific case, we cannot cluster b/c we only have data for a single year: 2015 (the only other two cases of a PJ with no copartisan presi at t+4 are 1983->1987 and 2019->2023, but in neither case do we have data on outcomes)
  }
  
  pj_copartisan_clustyear <- rd_hte (all, hi, lo, out_main[i], "PJ", "Copartisan", "Opposition"
                                     , mean_y_all = mean (filter (dpj, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                     , mean_y_hi = ifelse (
                                       i %in% grep ("_t2", out_main)
                                       , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                       , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                                     , mean_y_lo = ifelse (
                                       i %in% grep ("_t2", out_main)
                                       , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                       , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## UCR
  ducr$outcome_tmp <- unlist (ducr[,which (colnames (ducr) == out_main[i])])
  
  all <- with (filter (ducr, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=factor(year)))
  if (i %in% grep ("_t2", out_main)) { ## outcome measured at t+2
    hi <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=factor(year)));
    lo <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=factor(year)))
  } else {
    hi <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL)); ## opposite problem here
    lo <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=factor(year)))
  }
  
  ucr_copartisan_clustyear <- rd_hte (all, hi, lo, out_main[i], "UCR", "Copartisan", "Opposition"
                                      , mean_y_all = mean (filter (ducr, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                      , mean_y_hi = ifelse (
                                        i %in% grep ("_t2", out_main)
                                        , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                        , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                                      , mean_y_lo = ifelse (
                                        i %in% grep ("_t2", out_main)
                                        , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                        , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_copartisan_clustyear <- rbind (rd_results_copartisan_clustyear, cbind (matrix (unlist (
    list (pj_copartisan_clustyear, ucr_copartisan_clustyear)), nrow=8)))
}
rd_results_copartisan_clustyear


## favored copartisan president: PJ only
for (i in 1:length (out_main)){
  
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_main[i])])
  
  all <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=factor(year)))
  if (i %in% grep ("_t2", out_main)) { ## outcome measured at t+2
    hi <- with (filter (dpj, midterm=="concurrent" & favored_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=factor(year)));
    lo <- with (filter (dpj, midterm=="concurrent" & disfavored_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=factor(year)))
  } else {
    hi <- with (filter (dpj, midterm=="concurrent" & favored_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=factor(year)));
    lo <- with (filter (dpj, midterm=="concurrent" & disfavored_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=factor(year)))
  }
  
  pj_favored_clustyear <- rd_hte (all, hi, lo, out_main[i], "PJ", "Favored", "Disfavored"
                                  , mean_y_all = mean (filter (dpj, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                  , mean_y_hi = ifelse (
                                    i %in% grep ("_t2", out_main)
                                    , mean (filter (dpj, midterm=="concurrent" & favored_presi_t2==1 & !is.na (favored_presi_t2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                    , mean (filter (dpj, midterm=="concurrent" & favored_presi_t4==1 & !is.na (favored_presi_t4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                                  , mean_y_lo = ifelse (
                                    i %in% grep ("_t2", out_main)
                                    , mean (filter (dpj, midterm=="concurrent" & disfavored_presi_t2==0 & !is.na (disfavored_presi_t2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                    , mean (filter (dpj, midterm=="concurrent" & disfavored_presi_t4==0 & !is.na (disfavored_presi_t4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_favored_clustyear <- rbind (rd_results_favored_clustyear, cbind (matrix (unlist (
    list (pj_favored_clustyear)), nrow=8)))
}
rd_results_favored_clustyear


## (5.3.2) Table A7: full sample ####
tabFullSE <- rbind (
 matrix (rd_results_copartisan_clustyear[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),1]
         , byrow=F, nrow=6)
 , matrix (rd_results_copartisan_clustyear[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),5]
           , byrow=F, nrow=6))
(tabFullSE <- cbind (tabFullSE[,1:2], rep ("", nrow (tabFullSE)), tabFullSE[,3:4], rep ("", nrow (tabFullSE)), tabFullSE[,5:6], rep ("", nrow (tabFullSE)), tabFullSE[,7:8]))

rows <- rep ( c ("estimate ($\\hat{\\tau}_{\\textsc{rd}}$)"
         , "95\\% \\textsc{ci}"
         , "$p$-value"
         , "bwd."
         , "$N$"
         , "control mean"), 2)

Header1 <- paste ("\\toprule & \\multicolumn{2}{c}{} & & \\multicolumn{2}{c}{\\emph{vote share}} & & \\multicolumn{2}{c}{\\emph{seat share}} & & \\multicolumn{2}{c}{\\emph{vote share}} \\\\ \n")
Header2 <- paste ("& \\multicolumn{2}{c}{\\emph{winner}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(federal)}} \\\\ \\cmidrule{2-3} \\cmidrule{5-6} \\cmidrule{8-9} \\cmidrule{11-12} \n")
Header3 <- paste ("\\multicolumn{1}{l}{(a) \\textsc{pj}} & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{12}{l}{(b) \\textsc{ucr}} \\\\ \\midrule \n")

addtorow <- list()
addtorow$pos <- list()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 6
addtorow$command <- c (Header1, Header2, Header3, Header4)
print (xtable ( cbind (rows, tabFullSE)
        , align=c("l","l","c","c","c","c","c","c","c","c","c","c","c")
        , digits=2
        , caption="{Robustness checks (\\textsc{ii}): Clustering standard errors by year}"
        , label="T:robSE")
    , sanitize.text.function=function(x){x}
    , floating=TRUE
    , table.placement="t"
    , caption.placement="top" 
    , latex.environments="center"
    , size="footnotesize"
    , include.colnames=FALSE
    , include.rownames=FALSE
    , hline.after = c ()
    , add.to.row=addtorow )



### (5.4) 2nd order polynomials ####

## (5.4.1) running the models ####
rd_results_copartisan_poly2 <- NULL; rd_results_favored_poly2 <- NULL
for (i in 1:length (out_main)){
  
  ## PJ
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_main[i])])
  
  all <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio))
  if (i %in% grep ("_t2", out_main)) { ## outcome measured at t+2
    hi <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio))
  } else {
    hi <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio))
  }
  
  pj_copartisan_poly2 <- rd_hte (all, hi, lo, out_main[i], "PJ", "Copartisan", "Opposition"
                                 , mean_y_all = mean (filter (dpj, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                 , mean_y_hi = ifelse (
                                   i %in% grep ("_t2", out_main)
                                   , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                   , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                                 , mean_y_lo = ifelse (
                                   i %in% grep ("_t2", out_main)
                                   , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                   , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## UCR
  ducr$outcome_tmp <- unlist (ducr[,which (colnames (ducr) == out_main[i])])
  
  all <- with (filter (ducr, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio))
  if (i %in% grep ("_t2", out_main)) { ## outcome measured at t+2
    hi <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio));
    lo <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio))
  } else {
    hi <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio));
    lo <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio))
  }
  
  ucr_copartisan_poly2 <- rd_hte (all, hi, lo, out_main[i], "UCR", "Copartisan", "Opposition"
                                  , mean_y_all = mean (filter (ducr, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                  , mean_y_hi = ifelse (
                                    i %in% grep ("_t2", out_main)
                                    , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                    , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                                  , mean_y_lo = ifelse (
                                    i %in% grep ("_t2", out_main)
                                    , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                    , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_copartisan_poly2 <- rbind (rd_results_copartisan_poly2, cbind (matrix (unlist (
    list (pj_copartisan_poly2, ucr_copartisan_poly2)), nrow=8)))
}
rd_results_copartisan_poly2


## favored copartisan president: PJ only
for (i in 1:length (out_main)){
  
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_main[i])])
  
  all <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio))
  if (i %in% grep ("_t2", out_main)) { ## outcome measured at t+2
    hi <- with (filter (dpj, midterm=="concurrent" & favored_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & disfavored_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio))
  } else {
    hi <- with (filter (dpj, midterm=="concurrent" & favored_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & disfavored_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, p=2, q=3, cluster=municipio))
  }
  
  pj_favored_poly2 <- rd_hte (all, hi, lo, out_main[i], "PJ", "Favored", "Disfavored"
                              , mean_y_all = mean (filter (dpj, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                              , mean_y_hi = ifelse (
                                i %in% grep ("_t2", out_main)
                                , mean (filter (dpj, midterm=="concurrent" & favored_presi_t2==1 & !is.na (favored_presi_t2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                , mean (filter (dpj, midterm=="concurrent" & favored_presi_t4==1 & !is.na (favored_presi_t4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                              , mean_y_lo = ifelse (
                                i %in% grep ("_t2", out_main)
                                , mean (filter (dpj, midterm=="concurrent" & disfavored_presi_t2==0 & !is.na (disfavored_presi_t2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                , mean (filter (dpj, midterm=="concurrent" & disfavored_presi_t4==0 & !is.na (disfavored_presi_t4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_favored_poly2 <- rbind (rd_results_favored_poly2, cbind (matrix (unlist (
    list (pj_favored_poly2)), nrow=8)))
}
rd_results_favored_poly2


## (5.4.2) Table A8: full sample ####
tabFullPoly2 <- rbind (
 matrix (rd_results_copartisan_poly2[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),1]
         , byrow=F, nrow=6)
 , matrix (rd_results_copartisan_poly2[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),5]
           , byrow=F, nrow=6))
(tabFullPoly2 <- cbind (tabFullPoly2[,1:2], rep ("", nrow (tabFullPoly2)), tabFullPoly2[,3:4], rep ("", nrow (tabFullPoly2)), tabFullPoly2[,5:6], rep ("", nrow (tabFullPoly2)), tabFullPoly2[,7:8]))

rows <- rep ( c ("estimate ($\\hat{\\tau}_{\\textsc{rd}}$)"
         , "95\\% \\textsc{ci}"
         , "$p$-value"
         , "bwd."
         , "$N$"
         , "control mean"), 2)

Header1 <- paste ("\\toprule & \\multicolumn{2}{c}{} & & \\multicolumn{2}{c}{\\emph{vote share}} & & \\multicolumn{2}{c}{\\emph{seat share}} & & \\multicolumn{2}{c}{\\emph{vote share}} \\\\ \n")
Header2 <- paste ("& \\multicolumn{2}{c}{\\emph{winner}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(federal)}} \\\\ \\cmidrule{2-3} \\cmidrule{5-6} \\cmidrule{8-9} \\cmidrule{11-12} \n")
Header3 <- paste ("\\multicolumn{1}{l}{(a) \\textsc{pj}} & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{12}{l}{(b) \\textsc{ucr}} \\\\ \\midrule \n")

addtorow <- list()
addtorow$pos <- list()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 6
addtorow$command <- c (Header1, Header2, Header3, Header4)
print (xtable ( cbind (rows, tabFullPoly2)
        , align=c("l","l","c","c","c","c","c","c","c","c","c","c","c")
        , digits=2
        , caption="{Robustness checks (\\textsc{iii}): Employing second-order polynomials}"
        , label="T:robPoly2")
    , sanitize.text.function=function(x){x}
    , floating=TRUE
    , table.placement="t"
    , caption.placement="top" 
    , latex.environments="center"
    , size="footnotesize"
    , include.colnames=FALSE
    , include.rownames=FALSE
    , hline.after = c ()
    , add.to.row=addtorow )



### (5.5) Demeaned outcomes ####

## (5.5.1) running the models ####
rd_results_copartisan_fe <- NULL; rd_results_favored_fe <- NULL
for (i in 1:length (out_fe)){
  
  ## PJ
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_fe[i])])
  
  all <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_t2", out_fe)) { ## outcome measured at t+2
    hi <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  }
  
  pj_copartisan_fe <- rd_hte (all, hi, lo, out_fe[i], "PJ", "Copartisan", "Opposition"
                              , mean_y_all = mean (filter (dpj, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                              , mean_y_hi = ifelse (
                                i %in% grep ("_t2", out_fe)
                                , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                              , mean_y_lo = ifelse (
                                i %in% grep ("_t2", out_fe)
                                , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                , mean (filter (dpj, midterm=="concurrent" & copartisan_presi_t4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## UCR
  ducr$outcome_tmp <- unlist (ducr[,which (colnames (ducr) == out_fe[i])])
  
  all <- with (filter (ducr, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_t2", out_fe)) { ## outcome measured at t+2
    hi <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  }
  
  ucr_copartisan_fe <- rd_hte (all, hi, lo, out_fe[i], "UCR", "Copartisan", "Opposition"
                               , mean_y_all = mean (filter (ducr, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                               , mean_y_hi = ifelse (
                                 i %in% grep ("_t2", out_fe)
                                 , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                 , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                               , mean_y_lo = ifelse (
                                 i %in% grep ("_t2", out_fe)
                                 , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                 , mean (filter (ducr, midterm=="concurrent" & copartisan_presi_t4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_copartisan_fe <- rbind (rd_results_copartisan_fe, cbind (matrix (unlist (
    list (pj_copartisan_fe, ucr_copartisan_fe)), nrow=8)))
}
rd_results_copartisan_fe


## favored copartisan president: PJ only
for (i in 1:length (out_fe)){
  
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_fe[i])])
  
  all <- with (filter (dpj, midterm=="concurrent"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_t2", out_fe)) { ## outcome measured at t+2
    hi <- with (filter (dpj, midterm=="concurrent" & favored_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & disfavored_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (dpj, midterm=="concurrent" & favored_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="concurrent" & disfavored_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  }
  
  pj_favored_fe <- rd_hte (all, hi, lo, out_fe[i], "PJ", "Favored", "Disfavored"
                           , mean_y_all = mean (filter (dpj, midterm=="concurrent" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                           , mean_y_hi = ifelse (
                             i %in% grep ("_t2", out_fe)
                             , mean (filter (dpj, midterm=="concurrent" & favored_presi_t2==1 & !is.na (favored_presi_t2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                             , mean (filter (dpj, midterm=="concurrent" & favored_presi_t4==1 & !is.na (favored_presi_t4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                           , mean_y_lo = ifelse (
                             i %in% grep ("_t2", out_fe)
                             , mean (filter (dpj, midterm=="concurrent" & disfavored_presi_t2==0 & !is.na (disfavored_presi_t2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                             , mean (filter (dpj, midterm=="concurrent" & disfavored_presi_t4==0 & !is.na (disfavored_presi_t4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_favored_fe <- rbind (rd_results_favored_fe, cbind (matrix (unlist (
    list (pj_favored_fe)), nrow=8)))
}
rd_results_favored_fe


## (5.5.2) Table A9: full sample ####
tabFullFE <- rbind (
 matrix (rd_results_copartisan_fe[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),1]
         , byrow=F, nrow=6)
 , matrix (rd_results_copartisan_fe[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),5]
           , byrow=F, nrow=6))
(tabFullFE <- cbind (tabFullFE[,1:2], rep ("", nrow (tabFullFE)), tabFullFE[,3:4], rep ("", nrow (tabFullFE)), tabFullFE[,5:6], rep ("", nrow (tabFullFE)), tabFullFE[,7:8]))

rows <- rep ( c ("estimate ($\\hat{\\tau}_{\\textsc{rd}}$)"
         , "95\\% \\textsc{ci}"
         , "$p$-value"
         , "bwd."
         , "$N$"
         , "control mean"), 2)

Header1 <- paste ("\\toprule & \\multicolumn{2}{c}{} & & \\multicolumn{2}{c}{\\emph{vote share}} & & \\multicolumn{2}{c}{\\emph{seat share}} & & \\multicolumn{2}{c}{\\emph{vote share}} \\\\ \n")
Header2 <- paste ("& \\multicolumn{2}{c}{\\emph{winner}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(federal)}} \\\\ \\cmidrule{2-3} \\cmidrule{5-6} \\cmidrule{8-9} \\cmidrule{11-12} \n")
Header3 <- paste ("\\multicolumn{1}{l}{(a) \\textsc{pj}} & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{12}{l}{(b) \\textsc{ucr}} \\\\ \\midrule \n")

addtorow <- list()
addtorow$pos <- list()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 6
addtorow$command <- c (Header1, Header2, Header3, Header4)
print (xtable ( cbind (rows, tabFullFE)
        , align=c("l","l","c","c","c","c","c","c","c","c","c","c","c")
        , digits=2
        , caption="{Robustness checks (\\textsc{iv}): Demeaned outcomes}"
        , label="T:robFE")
    , sanitize.text.function=function(x){x}
    , floating=TRUE
    , table.placement="t"
    , caption.placement="top" 
    , latex.environments="center"
    , size="footnotesize"
    , include.colnames=FALSE
    , include.rownames=FALSE
    , hline.after = c ()
    , add.to.row=addtorow )



### (5.6) Placebos: effect of winning midterm election ####

## (5.6.1) running the models ####
rd_results_copartisan_midt <- NULL; rd_results_favored_midt <- NULL
for (i in 1:length (out_midt)){
  
  ## PJ
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_midt[i])])
  
  all <- with (filter (dpj, midterm=="midterm"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_t2", out_midt)) { ## outcome measured at t+2
    hi <- with (filter (dpj, midterm=="midterm" & copartisan_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="midterm" & copartisan_presi_t2==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (dpj, midterm=="midterm" & copartisan_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="midterm" & copartisan_presi_t4==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  }
  
  pj_copartisan_midt <- rd_hte (all, hi, lo, out_midt[i], "PJ", "Copartisan", "Opposition"
                                , mean_y_all = mean (filter (dpj, midterm=="midterm" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                , mean_y_hi = ifelse (
                                  i %in% grep ("_t2", out_midt)
                                  , mean (filter (dpj, midterm=="midterm" & copartisan_presi_t2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                  , mean (filter (dpj, midterm=="midterm" & copartisan_presi_t4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                                , mean_y_lo = ifelse (
                                  i %in% grep ("_t2", out_midt)
                                  , mean (filter (dpj, midterm=="midterm" & copartisan_presi_t2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                  , mean (filter (dpj, midterm=="midterm" & copartisan_presi_t4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## UCR
  ducr$outcome_tmp <- unlist (ducr[,which (colnames (ducr) == out_midt[i])])
  
  all <- with (filter (ducr, midterm=="midterm"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_t2", out_midt)) { ## outcome measured at t+2
    hi <- with (filter (ducr, midterm=="midterm" & copartisan_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (ducr, midterm=="midterm" & copartisan_presi_t2==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (ducr, midterm=="midterm" & copartisan_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (ducr, midterm=="midterm" & copartisan_presi_t4==0), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  }
  
  ucr_copartisan_midt <- rd_hte (all, hi, lo, out_midt[i], "UCR", "Copartisan", "Opposition"
                                 , mean_y_all = mean (filter (ducr, midterm=="midterm" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                 , mean_y_hi = ifelse (
                                   i %in% grep ("_t2", out_midt)
                                   , mean (filter (ducr, midterm=="midterm" & copartisan_presi_t2==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                   , mean (filter (ducr, midterm=="midterm" & copartisan_presi_t4==1 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                                 , mean_y_lo = ifelse (
                                   i %in% grep ("_t2", out_midt)
                                   , mean (filter (ducr, midterm=="midterm" & copartisan_presi_t2==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                                   , mean (filter (ducr, midterm=="midterm" & copartisan_presi_t4==0 & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_copartisan_midt <- rbind (rd_results_copartisan_midt, cbind (matrix (unlist (
    list (pj_copartisan_midt, ucr_copartisan_midt)), nrow=8)))
}
rd_results_copartisan_midt


## favored copartisan president: PJ only
for (i in 1:length (out_midt)){
  
  dpj$outcome_tmp <- unlist (dpj[,which (colnames (dpj) == out_midt[i])])
  
  all <- with (filter (dpj, midterm=="midterm"), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  if (i %in% grep ("_t2", out_midt)) { ## outcome measured at t+2
    hi <- with (filter (dpj, midterm=="midterm" & favored_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="midterm" & disfavored_presi_t2==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  } else {
    hi <- with (filter (dpj, midterm=="midterm" & favored_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio));
    lo <- with (filter (dpj, midterm=="midterm" & disfavored_presi_t4==1), rdrobust (y=outcome_tmp, x=running, bwselect="mserd", covs=NULL, cluster=municipio))
  }
  
  pj_favored_midt <- rd_hte (all, hi, lo, out_midt[i], "PJ", "Favored", "Disfavored"
                             , mean_y_all = mean (filter (dpj, midterm=="midterm" & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                             , mean_y_hi = ifelse (
                               i %in% grep ("_t2", out_midt)
                               , mean (filter (dpj, midterm=="midterm" & favored_presi_t2==1 & !is.na (favored_presi_t2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                               , mean (filter (dpj, midterm=="midterm" & favored_presi_t4==1 & !is.na (favored_presi_t4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp))
                             , mean_y_lo = ifelse (
                               i %in% grep ("_t2", out_midt)
                               , mean (filter (dpj, midterm=="midterm" & disfavored_presi_t2==0 & !is.na (disfavored_presi_t2) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)
                               , mean (filter (dpj, midterm=="midterm" & disfavored_presi_t4==0 & !is.na (disfavored_presi_t4) & !is.na (outcome_tmp) & !is.na (running) & running<0)$outcome_tmp)))
  
  
  ## exporting the results
  rd_results_favored_midt <- rbind (rd_results_favored_midt, cbind (matrix (unlist (
    list (pj_favored_midt)), nrow=8)))
}
rd_results_favored_midt


## (5.6.2) Table A5: full sample ####
tabFullMidt <- rbind (
 matrix (rd_results_copartisan_midt[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),1]
         , byrow=F, nrow=6)
 , matrix (rd_results_copartisan_midt[-c (1:2, 9:10, 17:18, 25:26, 33:34, 41:42, 49:50, 57:58),5]
           , byrow=F, nrow=6))
(tabFullMidt <- cbind (tabFullMidt[,1:2], rep ("", nrow (tabFullMidt)), tabFullMidt[,3:4], rep ("", nrow (tabFullMidt)), tabFullMidt[,5:6], rep ("", nrow (tabFullMidt)), tabFullMidt[,7:8]))

rows <- rep ( c ("estimate ($\\hat{\\tau}_{\\textsc{rd}}$)"
         , "95\\% \\textsc{ci}"
         , "$p$-value"
         , "bwd."
         , "$N$"
         , "control mean"), 2)

Header1 <- paste ("\\toprule & \\multicolumn{2}{c}{} & & \\multicolumn{2}{c}{\\emph{vote share}} & & \\multicolumn{2}{c}{\\emph{seat share}} & & \\multicolumn{2}{c}{\\emph{vote share}} \\\\ \n")
Header2 <- paste ("& \\multicolumn{2}{c}{\\emph{winner}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(municipal)}} & & \\multicolumn{2}{c}{\\emph{(federal)}} \\\\ \\cmidrule{2-3} \\cmidrule{5-6} \\cmidrule{8-9} \\cmidrule{11-12} \n")
Header3 <- paste ("\\multicolumn{1}{l}{(a) \\textsc{pj}} & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} & & \\multicolumn{1}{c}{$t+2$} & \\multicolumn{1}{c}{$t+4$} \\\\ \\midrule \n")
Header4 <- paste ("[2.5ex] \\multicolumn{12}{l}{(b) \\textsc{ucr}} \\\\ \\midrule \n")

addtorow <- list()
addtorow$pos <- list()
addtorow$pos[[1]] <- 0
addtorow$pos[[2]] <- 0
addtorow$pos[[3]] <- 0
addtorow$pos[[4]] <- 6
addtorow$command <- c (Header1, Header2, Header3, Header4)
print (xtable ( cbind (rows, tabFullMidt)
        , align=c("l","l","c","c","c","c","c","c","c","c","c","c","c")
        , digits=2
        , caption="{Robustness checks (\\textsc{v}): Placebo (midterm) elections}"
        , label="T:robMidt")
    , sanitize.text.function=function(x){x}
    , floating=TRUE
    , table.placement="t"
    , caption.placement="top" 
    , latex.environments="center"
    , size="footnotesize"
    , include.colnames=FALSE
    , include.rownames=FALSE
    , hline.after = c ()
    , add.to.row=addtorow )



### exporting everything ####
setwd (home)
save.image ("code/Results 2022-04-26.RData")

sink (NULL)
